{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Introduction to Quantitative Finance\n",
    "\n",
    "Copyright (c) 2019 Python Charmers Pty Ltd, Australia, <https://pythoncharmers.com>. All rights reserved.\n",
    "\n",
    "<img src=\"img/python_charmers_logo.png\" width=\"300\" alt=\"Python Charmers Logo\">\n",
    "\n",
    "Published under the Creative Commons Attribution-NonCommercial 4.0 International (CC BY-NC 4.0) license. See `LICENSE.md` for details.\n",
    "\n",
    "Sponsored by Tibra Global Services, <https://tibra.com>\n",
    "\n",
    "<img src=\"img/tibra_logo.png\" width=\"300\" alt=\"Tibra Logo\">\n",
    "\n",
    "\n",
    "## Module 1.5: Bayesian inference\n",
    "\n",
    "### 1.5.3 Representing prior knowledge: postcodes\n",
    "\n",
    "In the last module we looked at incorporating prior knowledge into our models, but effectively cheated by saying \"all outcomes are just as likely, until we have data\". Often we have more information than that.\n",
    "\n",
    "In this module we look at how to take prior information about the domain we are investigating, and using that information to alter our model.\n",
    "\n",
    "\n",
    "- France's *La Poste* has used automated sorting since 1964.\n",
    "- Handwritten digit recognition has been well studied.\n",
    "- In order to use digit recognition in practice for sorting mail, we need a *prior* model for how probable each postcode is, independent of each actual digitized hand-written digit image in front of us.\n",
    "\n",
    "A sample of some Australian and international postcodes:\n",
    "- 2000 (Sydney)\n",
    "- 3122 (Hawthorn, VIC)\n",
    "- 4350 (used for 44 towns near Toowoomba, QLD)\n",
    "- 8007 (PO boxes in Collins Street West)\n",
    "- A-1220 (Vienna, Austria)\n",
    "- Tsuen Wan (Hong Kong): no postcodes in HK\n",
    "- 02138 (Cambridge, MA)\n",
    "- EC1V 4AD (London)\n",
    "\n",
    "\n",
    "### A prior for Australian postcodes\n",
    "\n",
    "**What prior information do we have?**\n",
    "\n",
    "- Do all Australian postcodes have 4 digits? Yes.\n",
    "- What range? 0200 to 9944\n",
    "- States:\n",
    "   - NSW: postcodes 1000-1999 (PO boxes), 2000-2599, 2620-2899, 2921-2999\n",
    "   - ACT: 0200-0299 (PO boxes), 2600-2619, 2900-2920\n",
    "   - VIC: 3000-3999, 8000-8999 (PO boxes)\n",
    "   - QLD: 4000-4999, 9000-9999 (PO boxes)\n",
    "   - SA: 5000-5799, 5800-5999 (PO boxes)\n",
    "   - WA 6000-6797, 6800-6999 (PO boxes)\n",
    "   - TAS: 7000-7799, 7800-7999 (PO boxes)\n",
    "   - NT: 0800-0899, 0900-0999 (PO boxes)\n",
    "\n",
    "**Also:**\n",
    "\n",
    "- 25% of all mail goes to these CBD postcodes: 2000, 2001, 3000, 3001, 4000, 4001, 5000, 5001, 6000, 6001.\n",
    "- We can look up population for each postcode. Or, if we don't have population info by postcode, we could seed the prior with state population data.\n",
    "- Within each state xxxx, 80% of mail goes to x0xx and x1xx (metropolitan city areas and suburbs).\n",
    "\n",
    "### How do we encode this prior information for machine learning purposes?\n",
    "\n",
    "### Goal: construct a prior $p(\\textrm{postcode})$ over all 4-digit postcodes\n",
    "\n",
    "### Valid ranges\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "attributes": {
     "classes": [],
     "id": "",
     "n": "3"
    }
   },
   "outputs": [],
   "source": [
    "postcodes_by_state = dict((\n",
    "    ('Australian Capital Territory', set(range(2600, 2620)) | set(range(2900, 2920))),\n",
    "    ('New South Wales', set(range(2000, 3000)) - set(range(2600, 2620)) - set(range(2900, 2920))),\n",
    "    ('Victoria', set(range(3000, 4000))),\n",
    "    ('Queensland', set(range(4000, 5000))),\n",
    "    ('South Australia', set(range(5000, 5800))),\n",
    "    ('Western Australia', set(range(6000, 6798))),\n",
    "    ('Tasmania', set(range(7000, 7800))),\n",
    "    ('Northern Territory', set(range(800, 900)))\n",
    "))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### State populations\n",
    "\n",
    "We will start by using state populations as a proxy for really knowing the proportion of mail sent to each postcode.\n",
    "\n",
    "(If we obtain more data, we can update and improve our model by applying Bayes' theorem later.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "attributes": {
     "classes": [],
     "id": "",
     "n": "5"
    }
   },
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "\n",
    "state_populations = pd.read_hdf('data/aus_state_populations.h5')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "attributes": {
     "classes": [],
     "id": "",
     "n": "6"
    }
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Population</th>\n",
       "      <th>%</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>State</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>New South Wales</th>\n",
       "      <td>7757843</td>\n",
       "      <td>32.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Victoria</th>\n",
       "      <td>6100877</td>\n",
       "      <td>25.2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Queensland</th>\n",
       "      <td>4860448</td>\n",
       "      <td>20.1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>South Australia</th>\n",
       "      <td>1710804</td>\n",
       "      <td>7.1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Western Australia</th>\n",
       "      <td>2623164</td>\n",
       "      <td>10.8</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Tasmania</th>\n",
       "      <td>519783</td>\n",
       "      <td>2.1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Northern Territory</th>\n",
       "      <td>245657</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Australian Capital Territory</th>\n",
       "      <td>398349</td>\n",
       "      <td>1.6</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                              Population     %\n",
       "State                                         \n",
       "New South Wales                  7757843  32.0\n",
       "Victoria                         6100877  25.2\n",
       "Queensland                       4860448  20.1\n",
       "South Australia                  1710804   7.1\n",
       "Western Australia                2623164  10.8\n",
       "Tasmania                          519783   2.1\n",
       "Northern Territory                245657   1.0\n",
       "Australian Capital Territory      398349   1.6"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "state_populations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "These are the desired feature expectations for each state."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "attributes": {
     "classes": [],
     "id": "",
     "n": "7"
    }
   },
   "outputs": [],
   "source": [
    "# Source of the data:\n",
    "def fetch_state_populations():\n",
    "    url = 'http://www.ausstats.abs.gov.au/Ausstats/subscriber.nsf/0/D52DEAAFCEDF7B2ACA2580EB00133359/$File/31010do001_201609.xls'\n",
    "\n",
    "    state_pop = pd.read_excel(url, sheetname='Table_8', skiprows=6,\n",
    "                  names=['State', 'Population', '%'])\n",
    "\n",
    "    state_pop.set_index('State', inplace=True)\n",
    "\n",
    "    drop_row_idx = list(state_pop.index).index('Other Territories')\n",
    "\n",
    "    state_pop.drop(state_pop.index[drop_row_idx:], inplace=True)\n",
    "\n",
    "    state_pop['Population'] = state_pop['Population'].astype(int)\n",
    "    # state_pop.to_hdf('state_populations.h5', key='populations')\n",
    "    return state_pop"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### How to incorporate this?\n",
    "\n",
    "... to model the probability of e.g. $p(\\textrm{postcode}=3122)$?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "attributes": {
     "classes": [],
     "id": "",
     "n": "8"
    }
   },
   "outputs": [],
   "source": [
    "def prior_state(state):\n",
    "    return state_populations['%'].loc[state] / 100"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "attributes": {
     "classes": [],
     "id": "",
     "n": "9"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.32"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "prior_state('New South Wales')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we have a prior $p(\\text{state})$.\n",
    "\n",
    "### From the definition of conditional probability:\n",
    "\n",
    "$p(\\textrm{postcode}) = \\sum_{\\textrm{all states}} p(\\textrm{postcode | state}) p(\\textrm{state})$\n",
    "\n",
    "#### Exercise\n",
    "\n",
    "Assuming you have a function `prior_postcode_given_state(postcode, state)`, implement this as a function `prior_postcode(postcode)`.\n",
    "\n",
    "### Solution hint:\n",
    "\n",
    "Iterate over all states in `state_populations.index`.\n",
    "\n",
    "### Solution:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "attributes": {
     "classes": [],
     "id": "",
     "n": "10"
    }
   },
   "outputs": [],
   "source": [
    "def prior_postcode(postcode):\n",
    "    p = 0.0\n",
    "    for state in state_populations.index:\n",
    "        p += prior_postcode_given_state(postcode, state) * prior_state(state)\n",
    "    assert p <= 1\n",
    "    return p"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Exercise\n",
    "\n",
    "1. Write the function `prior_postcode_given_state(postcode, state)`.\n",
    "\n",
    "Assume you can assign equal probability to each valid postcode in the corresponding state (or 0 probability for the wrong state).\n",
    "\n",
    "\n",
    "You can test your code by trying these two examples:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.000252"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def prior_postcode_given_state(postcode, state):\n",
    "    if postcode in postcodes_by_state[state]:\n",
    "        return 1/len(postcodes_by_state[state])\n",
    "    else:\n",
    "        return 0\n",
    "\n",
    "\n",
    "\n",
    "prior_postcode_given_state(3122, 'Victoria')\n",
    "\n",
    "prior_postcode(3122)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*For solutions, see `solutions/prior_postcode_given_state.py`*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "%run -i solutions/prior_postcode_given_state.py"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### What did we do?\n",
    "\n",
    "We informally constructed a prior model that was as **flat** (uninformative) as possible **subject to a constraint** that the proportion of mail being delivered to a postcode is equal to the state's population, divided by the number of postcodes for that state.\n",
    "\n",
    "### Consider now: how would you update the model to reflect that ...\n",
    "\n",
    "1. 25% of all mail goes to one of the CBD postcodes; and\n",
    "2. Within each state xxxx, 80% of mail goes to x0xx and x1xx (metropolitan city areas and suburbs)?\n",
    "\n",
    "### Maximum entropy models: the easy way\n",
    "\n",
    "Here we see how to derive such prior models in a more systematic and principled way using the `maxentropy` package.\n",
    "\n",
    "You can install it like this:\n",
    "```\n",
    "pip install maxentropy\n",
    "```\n",
    "\n",
    "**Step 1: Set up the domain (or \"sample space\")**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "attributes": {
     "classes": [],
     "id": "",
     "n": "16"
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "samplespace = np.arange(10000, dtype=np.uint16)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "attributes": {
     "classes": [],
     "id": "",
     "n": "17"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([   0,    1,    2, ..., 9997, 9998, 9999], dtype=uint16)"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "samplespace"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Step 2: Set up a list of feature functions whose expectations you want to constrain**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "attributes": {
     "classes": [],
     "id": "",
     "n": "18"
    }
   },
   "outputs": [],
   "source": [
    "def is_valid(postcodes):\n",
    "    return [200 <= postcode < 10000 for postcode in postcodes]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "attributes": {
     "classes": [],
     "id": "",
     "n": "19"
    }
   },
   "outputs": [],
   "source": [
    "# def in_nsw(postcodes):\n",
    "#     return [postcode in postcodes_by_state['New South Wales'] for postcode in postcodes]\n",
    "# etc."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "attributes": {
     "classes": [],
     "id": "",
     "n": "20"
    }
   },
   "outputs": [],
   "source": [
    "def in_given_state(state):\n",
    "    def in_state(postcodes):\n",
    "        return [postcode in postcodes_by_state[state] for postcode in postcodes]\n",
    "    return in_state"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "attributes": {
     "classes": [],
     "id": "",
     "n": "21"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Index(['New South Wales', 'Victoria', 'Queensland', 'South Australia',\n",
       "       'Western Australia', 'Tasmania', 'Northern Territory',\n",
       "       'Australian Capital Territory'],\n",
       "      dtype='object', name='State')"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "state_populations.index"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "attributes": {
     "classes": [],
     "id": "",
     "n": "22"
    }
   },
   "outputs": [],
   "source": [
    "features = [is_valid] + \\\n",
    "           [in_given_state(state) for state in state_populations.index]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "attributes": {
     "classes": [],
     "id": "",
     "n": "23"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<function __main__.is_valid(postcodes)>,\n",
       " <function __main__.in_given_state.<locals>.in_state(postcodes)>,\n",
       " <function __main__.in_given_state.<locals>.in_state(postcodes)>,\n",
       " <function __main__.in_given_state.<locals>.in_state(postcodes)>,\n",
       " <function __main__.in_given_state.<locals>.in_state(postcodes)>,\n",
       " <function __main__.in_given_state.<locals>.in_state(postcodes)>,\n",
       " <function __main__.in_given_state.<locals>.in_state(postcodes)>,\n",
       " <function __main__.in_given_state.<locals>.in_state(postcodes)>,\n",
       " <function __main__.in_given_state.<locals>.in_state(postcodes)>]"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "features"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Step 3: create a `MinDivergenceModel` object from this list of features and sample space**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "attributes": {
     "classes": [],
     "id": "",
     "n": "25"
    }
   },
   "outputs": [
    {
     "ename": "ModuleNotFoundError",
     "evalue": "No module named 'maxentropy'",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mModuleNotFoundError\u001b[0m                       Traceback (most recent call last)",
      "Cell \u001b[1;32mIn[21], line 1\u001b[0m\n\u001b[1;32m----> 1\u001b[0m \u001b[38;5;28;01mfrom\u001b[39;00m \u001b[38;5;21;01mmaxentropy\u001b[39;00m \u001b[38;5;28;01mimport\u001b[39;00m MinDivergenceModel\n\u001b[0;32m      3\u001b[0m model \u001b[38;5;241m=\u001b[39m MinDivergenceModel(features, samplespace, vectorized\u001b[38;5;241m=\u001b[39m\u001b[38;5;28;01mTrue\u001b[39;00m)\n",
      "\u001b[1;31mModuleNotFoundError\u001b[0m: No module named 'maxentropy'"
     ]
    }
   ],
   "source": [
    "from maxentropy import MinDivergenceModel\n",
    "\n",
    "model = MinDivergenceModel(features, samplespace, vectorized=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Step 4: define your desired array of expected feature function values (one for each feature)**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "attributes": {
     "classes": [],
     "id": "",
     "n": "26"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "State\n",
       "New South Wales                 0.320\n",
       "Victoria                        0.252\n",
       "Queensland                      0.201\n",
       "South Australia                 0.071\n",
       "Western Australia               0.108\n",
       "Tasmania                        0.021\n",
       "Northern Territory              0.010\n",
       "Australian Capital Territory    0.016\n",
       "Name: %, dtype: float64"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pop = state_populations['%'] / 100\n",
    "pop"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "attributes": {
     "classes": [],
     "id": "",
     "n": "27"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "99.9"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "state_populations['%'].sum()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(This excludes the other territories, like Norfolk Island. Ignore this for now.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "attributes": {
     "classes": [],
     "id": "",
     "n": "28"
    }
   },
   "outputs": [],
   "source": [
    "# Target expectations\n",
    "k = np.atleast_2d(np.r_[1, pop.values])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "attributes": {
     "classes": [],
     "id": "",
     "n": "29"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[1.   , 0.32 , 0.252, 0.201, 0.071, 0.108, 0.021, 0.01 , 0.016]])"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "k"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "attributes": {
     "classes": [],
     "id": "",
     "n": "30"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(features) == k.shape[1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Step 5: fit your model under those constraints**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "attributes": {
     "classes": [],
     "id": "",
     "n": "31"
    }
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/ashkfarjami/opt/anaconda3/lib/python3.8/site-packages/sklearn/base.py:209: FutureWarning: From version 0.24, get_params will raise an AttributeError if a parameter cannot be retrieved as an instance attribute. Previously it would return None.\n",
      "  warnings.warn('From version 0.24, get_params will raise an '\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "MinDivergenceModel(features=[<function is_valid at 0x7fe5b3444e50>,\n",
       "                             <function in_given_state.<locals>.in_state at 0x7fe5b34444c0>,\n",
       "                             <function in_given_state.<locals>.in_state at 0x7fe5b345a040>,\n",
       "                             <function in_given_state.<locals>.in_state at 0x7fe5b345a0d0>,\n",
       "                             <function in_given_state.<locals>.in_state at 0x7fe5b345a160>,\n",
       "                             <function in_given_state.<locals>.in_state at 0x7fe5b345a1f0>,\n",
       "                             <function in_given_state.<locals>.in_state at 0x7fe5b345a280>,\n",
       "                             <function in_given_state.<locals>.in_state at 0x7fe5b345a310>,\n",
       "                             <function in_given_state.<locals>.in_state at 0x7fe5b345a3a0>],\n",
       "                   samplespace=array([   0,    1,    2, ..., 9997, 9998, 9999], dtype=uint16),\n",
       "                   vectorized=True)"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model.fit(k)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "attributes": {
     "classes": [],
     "id": "",
     "n": "32"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[-4.04218898e-08, -7.21599679e-08, -7.89633941e-07,\n",
       "         7.45644801e-07, -3.70129965e-07,  1.00904637e-06,\n",
       "         5.08948233e-08, -2.96428910e-07, -1.79097180e-07]])"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model.expectations() - k"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "attributes": {
     "classes": [],
     "id": "",
     "n": "33"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.allclose(model.expectations(), k, atol=1e-6)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Result: our fitted prior model is given by `model.probdist()`**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {
    "attributes": {
     "classes": [],
     "id": "",
     "n": "34"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([2.02108897e-10, 2.02108897e-10, 2.02108897e-10, ...,\n",
       "       2.32417815e-07, 2.32417815e-07, 2.32417815e-07])"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model.probdist()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {
    "attributes": {
     "classes": [],
     "id": "",
     "n": "36"
    }
   },
   "outputs": [],
   "source": [
    "assert len(model.probdist() == len(samplespace))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now have a prior probability $\\textrm{prior}(\\textrm{postcode})$ for each 4-digit postcode.\n",
    "\n",
    "### What are the most probable postcodes?\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {
    "attributes": {
     "classes": [],
     "id": "",
     "n": "37"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([2915, 2910, 2917, ...,  128,  127,    0])"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "p = model.probdist()\n",
    "np.argsort(p)[::-1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Visualized\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {
    "attributes": {
     "classes": [],
     "id": "",
     "n": "38"
    }
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {
    "attributes": {
     "classes": [],
     "id": "",
     "n": "39"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, 'probability')"
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAuoAAAE9CAYAAABdtEesAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAm6klEQVR4nO3df7xddX3n+9ebRIRqMSCBwQQKllTLjzuIR0I7aG25lAS0iW2ZQagwDPOgWDOtM9NbIw+5Xjs+eKS3nepQKRQVTdqhKWOLZBSGQUb06gjmpCAQKHJEC4ekEO5cEUGlkc/9Y69jdw8nZ+8kZ2evc87r+Xisx9rru77f7/ouljnn7TrftXaqCkmSJEntst+wByBJkiTpxQzqkiRJUgsZ1CVJkqQWMqhLkiRJLWRQlyRJklrIoC5JkiS10MJhD6CtDj300Dr66KOHPQxJkiTNYVu2bHmqqhZPtc+gvgtHH300o6Ojwx6GJEmS5rAkf7urfU59kSRJklrIoC5JkiS1kEFdkiRJaiGDuiRJktRCBnVJkiSphQzqkiRJUgsZ1CVJkqQWGmhQT7IiyUNJxpKsnWJ/klzZ7L83ycm92iY5JMltSR5u1gdP6vOoJN9N8ttdZa9Pcl/T15VJMqhzliRJkmbCwL7wKMkC4CrgDGAc2JxkU1U90FVtJbCsWZYDVwPLe7RdC9xeVeuaAL8WeE9Xnx8Cbpk0nKuBS4A7gZuBFVPUkeal6+96lMtuvO9H21e87UTOW37UHvV17GWfZecLnc8nLX0Fn15z2kwMUZKkeWmQd9RPAcaq6pGqeh7YCKyaVGcVsKE67gQWJTmiR9tVwPrm83pg9URnSVYDjwBbu8qOAA6qqq9UVQEbuttI89nkkA5w2Y33cf1dj+52X90hHeCe8adZ/ZEv7e0QJUmatwYZ1JcAj3Vtjzdl/dSZru3hVbUdoFkfBpDkZXTurH9gimOM9xgHTR+XJBlNMrpjx45pT06aC265f/tulU+nO6RPuH/bd3a7H0mS1DHIoD7VPPDqs04/bSf7APChqvruHoyjU1h1bVWNVNXI4sWLexxOmv1WnnDEbpVPZ+EUP01OeNVBu92PJEnqGGRQHweO7NpeCmzrs850bZ9oprNMTGt5silfDvzfSb4FvBu4LMmapq+lPcYhzUvnLT+KK9524j8q29M56mNXnP2Pwrpz1CVJ2jsDe5gU2AwsS3IM8DhwLnDepDqbgDVJNtIJ2k9X1fYkO6Zpuwm4EFjXrG8CqKo3TnSa5P8CvltVH2m2n0lyKnAXcAHwRzN/utLsdN7yo/b44dHJxq44e0b6kSRJAwzqVbWzuaN9K7AAuK6qtia5tNl/DZ03sJwFjAHPARdN17bpeh1wQ5KLgUeBc/oYzjuBTwIH0nnbi298kSRJUqul8yIUTTYyMlKjo6PDHoYkSZLmsCRbqmpkqn1+M6kkSZLUQgZ1SZIkqYUM6pIkSVILGdQlSZKkFjKoS5IkSS1kUJckSZJayKAuSZIktZBBXZIkSWohg7okSZLUQgZ1SZIkqYUM6pIkSVILGdQlSZKkFjKoS5IkSS1kUJckSZJayKAuSZIktZBBXZIkSWohg7okSZLUQgZ1SZIkqYUM6pIkSVILGdQlSZKkFjKoS5IkSS000KCeZEWSh5KMJVk7xf4kubLZf2+Sk3u1TXJIktuSPNysD27KT0lyT7N8Lcnbutrc0fQ1sf+wQZ63JEmStLcGFtSTLACuAlYCxwFvT3LcpGorgWXNcglwdR9t1wK3V9Uy4PZmG+B+YKSqTgJWAH+SZGHXsc6vqpOa5ckZPVlJkiRphg3yjvopwFhVPVJVzwMbgVWT6qwCNlTHncCiJEf0aLsKWN98Xg+sBqiq56pqZ1N+AFADOi9JkiRp4AYZ1JcAj3Vtjzdl/dSZru3hVbUdoFn/aBpLkuVJtgL3AZd2BXeATzTTXi5Pkj0/LUmSJGnwBhnUpwrDk+9y76pOP21fXKHqrqo6HngD8N4kBzS7zq+qE4E3Nss7phxwckmS0SSjO3bs6HU4SZIkaWAGGdTHgSO7tpcC2/qsM13bJ5rpMTTrF803r6oHgWeBE5rtx5v1M8D1dKbWvEhVXVtVI1U1snjx4j5OUZIkSRqMQQb1zcCyJMck2R84F9g0qc4m4ILm7S+nAk8301mma7sJuLD5fCFwE0BTd2Hz+SeA1wDfSrIwyaFN+UuAt9B58FSSJElqrYW9q+yZqtqZZA1wK7AAuK6qtia5tNl/DXAzcBYwBjwHXDRd26brdcANSS4GHgXOacpPA9Ym+XvgBeA3quqpJC8Dbm1C+gLgc8BHB3XekiRJ0kxIlS9HmcrIyEiNjo4OexiSJEmaw5JsqaqRqfb5zaSSJElSCxnUJUmSpBYyqEuSJEktZFCXJEmSWmhgb32RNLWTPnAr3/7ezt4V56DVJ72KD5/7umEPQ5KkWcE76tI+NJ9DOsCn79nGuzfePexhSJI0KxjUpX1oPof0CXd8fcewhyBJ0qxgUJf2oUUHOtvszT+1eNhDkCRpVjCoS/vQPe8/c16HdeeoS5LUv/mbGKQhuef9Zw57CJIkaRbwjrokSZLUQgZ1SZIkqYUM6pIkSVILGdQlSZKkFjKoS5IkSS1kUJckSZJayKAuSZIktZBBXZIkSWohg7okSZLUQgZ1SZIkqYUM6pIkSVILGdQlSZKkFhpoUE+yIslDScaSrJ1if5Jc2ey/N8nJvdomOSTJbUkebtYHN+WnJLmnWb6W5G1dbV6f5L6mryuTZJDnLUmSJO2tgQX1JAuAq4CVwHHA25McN6naSmBZs1wCXN1H27XA7VW1DLi92Qa4HxipqpOAFcCfJFnY7Lu66X/iWCtm9GQlSZKkGTbIO+qnAGNV9UhVPQ9sBFZNqrMK2FAddwKLkhzRo+0qYH3zeT2wGqCqnquqnU35AUABNP0dVFVfqaoCNky0kSRJktpqkEF9CfBY1/Z4U9ZPnenaHl5V2wGa9WETlZIsT7IVuA+4tAnuS5r2041DkiRJapVBBvWp5oFXn3X6afviClV3VdXxwBuA9yY5YHf6SnJJktEkozt27Oh1OEmSJGlgBhnUx4Eju7aXAtv6rDNd2yea6SwT01qenHzgqnoQeBY4oelraY9xTLS7tqpGqmpk8eLF056cJEmSNEiDDOqbgWVJjkmyP3AusGlSnU3ABc3bX04Fnm6ms0zXdhNwYfP5QuAmgKbuwubzTwCvAb7V9PdMklObt71cMNFGkiRJaquFvavsmaramWQNcCuwALiuqrYmubTZfw1wM3AWMAY8B1w0Xdum63XADUkuBh4FzmnKTwPWJvl74AXgN6rqqWbfO4FPAgcCtzSLJEmS1FrpvAhFk42MjNTo6OiwhyFJkqQ5LMmWqhqZap/fTCpJkiS1kEFdkiRJaiGDuiRJktRCBnVJkiSphQzqkiRJUgsZ1CVJkqQWMqhLkiRJLWRQlyRJklrIoC5JkiS1kEFdkiRJaiGDuiRJktRCBnVJkiSphQzqkiRJUgsZ1CVJkqQWMqhLkiRJLWRQlyRJklrIoC5JkiS1kEFdkiRJaiGDuiRJktRCBnVJkiSphQzqkiRJUgsZ1CVJkqQWGmhQT7IiyUNJxpKsnWJ/klzZ7L83ycm92iY5JMltSR5u1gc35Wck2ZLkvmb9C11t7mj6uqdZDhvkeUuSJEl7a+GgOk6yALgKOAMYBzYn2VRVD3RVWwksa5blwNXA8h5t1wK3V9W6JsCvBd4DPAW8taq2JTkBuBVY0nWs86tqdFDnK2nXfvp9t/C9nS8MexhDccXbTuS85UcNexiSpFlokHfUTwHGquqRqnoe2AismlRnFbChOu4EFiU5okfbVcD65vN6YDVAVd1dVdua8q3AAUleOqBzk9Sn+RzSAS678T6uv+vRYQ9DkjQLDTKoLwEe69oe5x/f4Z6uznRtD6+q7QDNeqppLL8C3F1VP+gq+0Qz7eXyJNndk5G0Z+ZzSJ9wy/3bhz0ESdIsNMigPlUYrj7r9NN26oMmxwO/B/x6V/H5VXUi8MZmeccu2l6SZDTJ6I4dO/o5nKQeDlzoM+srTzhi2EOQJM1Cg/wNOg4c2bW9FNjWZ53p2j7RTI+hWT85USnJUuBG4IKq+sZEeVU93qyfAa6nM7XmRarq2qoaqaqRxYsX93makqbz4AdXzuuw7hx1SdKeGtjDpMBmYFmSY4DHgXOB8ybV2QSsSbKRzsOkT1fV9iQ7pmm7CbgQWNesbwJIsgj4LPDeqvryxAGSLAQWVdVTSV4CvAX43ADOV9IuPPjBlcMegiRJs87AgnpV7Uyyhs7bVxYA11XV1iSXNvuvAW4GzgLGgOeAi6Zr23S9DrghycXAo8A5Tfka4Fjg8iSXN2W/CDwL3NqE9AV0QvpHB3XekiRJ0kxIVV9Tv+edkZGRGh31bY6SJEkanCRbqmpkqn3zd+KoJEmS1GIGdUmSJKmF+grqSf4yydlJDPaSJEnSPtBv8L6azltXHk6yLslrBzgmSZIkad7rK6hX1eeq6nzgZOBbwG1J/meSi5q3qUiSJEmaQX1PZUnySuBfAv8auBv4T3SC+20DGZkkSZI0j/X1HvUkfwW8FvhT4K1Vtb3Z9RdJfIehJEmSNMP6/cKjj1XVzd0FSV5aVT/Y1XsfJUmSJO25fqe+fHCKsq/M5EAkSZIk/YNp76gn+SfAEuDAJK8D0uw6CPixAY9NkiRJmrd6TX05k84DpEuBP+wqfwa4bEBjkiRJkua9aYN6Va0H1if5lar6y300JkmSJGne6zX15deq6s+Ao5P8u8n7q+oPp2gmSZIkaS/1mvrysmb98kEPRJIkSdI/6DX15U+a9Qf2zXAkSZIkQe+pL1dOt7+qfnNmhyNJkiQJek992bJPRiFJkiTpH+nnrS+SJEmS9rFeU18+XFXvTvJfgZq8v6p+aWAjkyRJkuaxXlNf/rRZ/8GgByJJkiTpH/Sa+rKlWX8hyf7Aa+ncWX+oqp7fB+OTJEmS5qVed9QBSHI2cA3wDSDAMUl+vapuGeTgJEmSpPlqvz7r/Ufg56vqzVX1c8DPAx/q1SjJiiQPJRlLsnaK/UlyZbP/3iQn92qb5JAktyV5uFkf3JSfkWRLkvua9S90tXl9Uz7WHC99nrckSZI0FH3dUQeerKqxru1HgCena5BkAXAVcAYwDmxOsqmqHuiqthJY1izLgauB5T3argVur6p1TYBfC7wHeAp4a1VtS3ICcCuwpDnO1cAlwJ3AzcAKwL8GSBqYYy/7LDtfGPYo9r39gEfWnT3sYUjSnDDtHfUkv5zkl4GtSW5O8i+TXAj8V2Bzj75PAcaq6pFmPvtGYNWkOquADdVxJ7AoyRE92q4CJl4buR5YDVBVd1fVtqZ8K3BAkpc2/R1UVV+pqgI2TLSRpEGYryEd4AXg1Ws/O+xhSNKc0OuO+lu7Pj8B/FzzeQdwcI+2S4DHurbH6dw171VnSY+2h1fVdoCq2p7ksCmO/SvA3VX1gyRLmvaTjyFJAzFfQ/qEeX76kjRjer315aK96HuqeeCT38W+qzr9tJ36oMnxwO8Bv7gb45hoewmdKTIcddRR/RxOkl5k4X7zO6z3+/CTJGl6ff08TXJAkncl+eMk100sPZqNA0d2bS8FtvVZZ7q2TzTTWWjWP5orn2QpcCNwQVV9o+sYS3uMA4CquraqRqpqZPHixT1OT5KmNnbF2Sycp2nVOeqSNHP6fZj0T4G/Ac4Efhc4H3iwR5vNwLIkxwCPA+cC502qswlYk2QjnaktTzfTWXZM03YTcCGwrlnfBJBkEfBZ4L1V9eWJAzT9PZPkVOAu4ALgj/o8b0naI2NXGFYlSXun33s+x1bV5cCzVbUeOBs4cboGVbUTWEPn7SsPAjdU1dYklya5tKl2M503yIwBHwV+Y7q2TZt1wBlJHqbzVph1Tfka4Fjg8iT3NMvE/PV3Ah9rjvMNfOOLJEmSWi6dF6H0qJR8tapOSfJFOmH674CvVtWrBz3AYRkZGanR0dFhD0OSJElzWJItVTUy1b5+p75c23yx0OV0pp68vPksSZIkaQD6CupV9bHm4xeAOXsXXZIkSWqLft/68sokf5Tkr5NsSfLhJK8c9OAkSZKk+arfh0k30nkN4q8Avwo8BfzFoAYlSZIkzXf9zlE/pKr+Q9f2B5OsHsB4JEmSJNH/HfXPJzk3yX7N8s/pvLNckiRJ0gBMe0c9yTNAAQH+HfBnza79gO8C7x/o6CRJkqR5atqgXlU/vq8GIkmSJOkf9DtHnSS/BLyp2byjqj4zmCFJkiRJ6vf1jOuA3wIeaJbfasokSZIkDUC/d9TPAk6qqhcAkqwH7gbWDmpgkiRJ0nzW71tfABZ1fX7FDI9DkiRJUpd+76hfAdyd5PN03gDzJuC9AxuVJEmSNM/1DOpJ9gNeAE4F3kAnqL+nqv5uwGOTJEmS5q2eQb2qXkiypqpuADbtgzFJkiRJ816/c9RvS/LbSY5McsjEMtCRSZIkSfNYv3PU/xWdbyj9jUnlr57Z4UiSJEmC/oP6cXRC+ml0Avv/A1wzqEFJkiRJ812/QX098B3gymb77U3ZPx/EoCRJkqT5rt+g/pqq+qdd259P8rVBDEiSJElS/w+T3p3k1ImNJMuBLw9mSJIkSZL6vaO+HLggyaPN9lHAg0nuA6qq/reBjE6SJEmap/q9o74COAb4uWY5BjgLeAvw1l01SrIiyUNJxpKsnWJ/klzZ7L83ycm92javhrwtycPN+uCm/JVJPp/ku0k+Muk4dzR93dMsh/V53pIkSdJQ9HVHvar+dnc7TrIAuAo4AxgHNifZVFUPdFVbCSxrluXA1cDyHm3XArdX1bomwK8F3gN8H7gcOKFZJju/qkZ39zwkSZKkYej3jvqeOAUYq6pHqup5YCOwalKdVcCG6rgTWJTkiB5tV9F54wzNejVAVT1bVV+iE9glSZKkWW2QQX0J8FjX9nhT1k+d6doeXlXbAZp1v9NYPtFMe7k8SfpsI0mSJA3FIIP6VGG4+qzTT9vdcX5VnQi8sVneMVWlJJckGU0yumPHjr04nCRJkrR3BhnUx4Eju7aXAtv6rDNd2yea6TE06yd7DaSqHm/WzwDX05laM1W9a6tqpKpGFi9e3KtbSZIkaWAGGdQ3A8uSHJNkf+BcYNOkOpvovPYxzXvan26ms0zXdhNwYfP5QuCm6QaRZGGSQ5vPL6Hzppr79/70JEmSpMHp9z3qu62qdiZZA9wKLACuq6qtSS5t9l8D3EznNY9jwHPARdO1bbpeB9yQ5GLgUeCciWMm+RZwELB/ktXALwJ/C9zahPQFwOeAjw7qvCVJkqSZkKq9mfo9d42MjNToqG9zlCRJ0uAk2VJVI1PtG+TUF0mSJEl7yKAuSZIktdDA5qhLkjSXXX/Xo1x2433DHsaMWLroAL609vRhD0PSJN5RlyRpN82lkA4w/u3vc9q624c9DEmTGNQlSdpNt9y/fdhDmHGPf/v7wx6CpEkM6pIk7aaVJxwx7CHMuCWLDhj2ECRNYlCXJGk3nbf8KK5424nDHsaMcY661E4+TCpJ0h44b/lRnLf8qGEPQ9Ic5h11SZIkqYUM6pIkSVILGdQlSZKkFjKoS5IkSS1kUJckSZJayKAuSZIktZBBXZIkSWohg7okSZLUQgZ1SZIkqYUM6pIkSVILGdQlSZKkFjKoS5IkSS1kUJckSZJayKAuSZIktdBAg3qSFUkeSjKWZO0U+5Pkymb/vUlO7tU2ySFJbkvycLM+uCl/ZZLPJ/luko9MOs7rk9zX9HVlkgzyvCVJkqS9NbCgnmQBcBWwEjgOeHuS4yZVWwksa5ZLgKv7aLsWuL2qlgG3N9sA3wcuB357iuFc3fQ/cawVM3CKkiRJ0sAsHGDfpwBjVfUIQJKNwCrgga46q4ANVVXAnUkWJTkCOHqatquANzft1wN3AO+pqmeBLyU5tnsQTX8HVdVXmu0NwGrglhk+XwEn/J//je8+/8Nd7l+2+GXc9u/fvO8GJEmSNEsNcurLEuCxru3xpqyfOtO1PbyqtgM068P6GMd4j3FoBvQK6QAP73iWM/7jHftmQJIkSbPYIIP6VPPAq886/bSdyXF0KiaXJBlNMrpjx449PNz81SukT/jGU88OeCSSJEmz3yCD+jhwZNf2UmBbn3Wma/tEM51lYlrLk32MY2mPcQBQVddW1UhVjSxevLhHt5rs5fsv6KveTx76sgGPRJIkafYbZFDfDCxLckyS/YFzgU2T6mwCLmje/nIq8HQznWW6tpuAC5vPFwI3TTeIpr9nkpzavO3lgl5ttGfu/90VPcO6c9QlSZL6M7CHSatqZ5I1wK3AAuC6qtqa5NJm/zXAzcBZwBjwHHDRdG2brtcBNyS5GHgUOGfimEm+BRwE7J9kNfCLVfUA8E7gk8CBdB4i9UHSAbn/d32hjiRJ0kxI54UrmmxkZKRGR0eHPQxJkiTNYUm2VNXIVPsG+XpGSdI88uq1n+WFYQ9iBhy4cD8e/ODKYQ9Dkgb7zaSSpPlhroR0gO/tfIGffp8zJCUNn0FdkrTX5kpIn/C9nXPtjCTNRgZ1SdJem2u/TA5cONfOSNJs5E8iSdJee2Td2XPmF4pz1CW1hQ+TSpJmxCPrzh72ECRpTpkrN0AkSZKkOcWgLkmSJLWQQV2SJElqIYO6JEmS1EIGdUmSJKmFDOqSJElSCxnUJUmSpBYyqEuSJEktZFCXJEmSWsigLkmSJLWQQV2SJElqIYO6JEmS1EIGdUmSJKmFDOqSJElSCxnUJUmSpBYyqEuSJEktNNCgnmRFkoeSjCVZO8X+JLmy2X9vkpN7tU1ySJLbkjzcrA/u2vfepv5DSc7sKr+jKbunWQ4b5HlLkiRJe2tgQT3JAuAqYCVwHPD2JMdNqrYSWNYslwBX99F2LXB7VS0Dbm+2afafCxwPrAD+uOlnwvlVdVKzPDnT5ytJkiTNpEHeUT8FGKuqR6rqeWAjsGpSnVXAhuq4E1iU5IgebVcB65vP64HVXeUbq+oHVfVNYKzpR5IkSZp1BhnUlwCPdW2PN2X91Jmu7eFVtR2gWU9MY+l1vE80014uT5LdPx1JkiRp3xlkUJ8qDFefdfppuzvHO7+qTgTe2CzvmLKD5JIko0lGd+zY0eNwkiRJ0uAMMqiPA0d2bS8FtvVZZ7q2TzTTY2jWE/PNd9mmqh5v1s8A17OLKTFVdW1VjVTVyOLFi/s4RUmSJGkwBhnUNwPLkhyTZH86D3pumlRnE3BB8/aXU4Gnm+ks07XdBFzYfL4QuKmr/NwkL01yDJ0HVL+aZGGSQwGSvAR4C3D/IE5YkiRJmikLB9VxVe1Msga4FVgAXFdVW5Nc2uy/BrgZOIvOg5/PARdN17bpeh1wQ5KLgUeBc5o2W5PcADwA7ATeVVU/TPIy4NYmpC8APgd8dFDnLUmSJM2EVPWa+j0/jYyM1Ojo6LCHIUmSpDksyZaqGplqn99MKkmSJLWQQV2SJElqIYO6JEmS1EIGdUmSJKmFDOqSJElSCxnUJUmSpBYyqEuSJEktNLAvPJIkSXPLuzfezafv2TbsYcyIpYsO4EtrTx/2MKRpeUddkiT1NJdCOsD4t7/PaetuH/YwpGkZ1CVJUk93fH3HsIcw4x7/9veHPQRpWgZ1SZLU05t/avGwhzDjliw6YNhDkKZlUJckST19+NzXsfqkVw17GDPGOeqaDVJVwx5DK42MjNTo6OiwhyFJkqQ5LMmWqhqZap931CVJkqQWMqhLkiRJLWRQlyRJklrIoC5JkiS1kEFdkiRJaiGDuiRJktRCBnVJkiSphQzqkiRJUgsZ1CVJkqQWGmhQT7IiyUNJxpKsnWJ/klzZ7L83ycm92iY5JMltSR5u1gd37XtvU/+hJGd2lb8+yX3NviuTZJDnLUmSJO2thYPqOMkC4CrgDGAc2JxkU1U90FVtJbCsWZYDVwPLe7RdC9xeVeuaAL8WeE+S44BzgeOBVwGfS/JTVfXDpt9LgDuBm4EVwC2DOvc9tfojX+Ke8aeHPYwZ8aZlh7Lh4uU/2j7pA7fy7e/tHOKI2uekpa/g02tOG/YwJDUu+PhdfPHhp3a73cL9YOyKswcwosF7wwdvY8d3n9/tdgcu3I8HP7hyACMavD09Z80PbfvdPMg76qcAY1X1SFU9D2wEVk2qswrYUB13AouSHNGj7SpgffN5PbC6q3xjVf2gqr4JjAGnNP0dVFVfqaoCNnS1aY25FNIBvvjwU1zw8bsAQ/qu3DP+NKs/8qVhD0MSex7SAXa+AMde9tkZHtHg7U1g/d7OF/jp97XufldPhnT10rbfzYMM6kuAx7q2x5uyfupM1/bwqtoO0KwP66Ov8R7jACDJJUlGk4zu2LFj2pObafdv+84+Pd6+8NVv/S8AQ/o05uJ1l2ajiZ9Xe2rnCzM0kH1obwPr92bhSRvS1Y82/W4eZFCfah549Vmnn7b9Hq/vvqrq2qoaqaqRxYsX9zjczDrhVQft0+PtC6ccfQgAiw4c2AyrWW8uXndpNpr4ebWnFs7CVzMsfvn+e9X+wFl40nt7zpof2vS7eZD/ysaBI7u2lwLb+qwzXdsnmuksNOsn++hraY9xDN2n15zGSUtfMexhzJjuOer3vP9Mw/oU2jYPTprPNly8nDctO3SP2s7WOeqb33fGHgfX2TpHfW/OWfND2343pzNtewAdJwuBrwOnA48Dm4HzqmprV52zgTXAWXQeJr2yqk6Zrm2S3wf+366HSQ+pqt9JcjxwPZ357a8CbgeWVdUPk2wG/g1wF52HSf+oqm6ebvwjIyM1Ojo6Y/89JEmSpMmSbKmqkan2Dew2Z1XtTLIGuBVYAFzXBO1Lm/3X0AnNZ9F58PM54KLp2jZdrwNuSHIx8ChwTtNma5IbgAeAncC7mje+ALwT+CRwIJ23vcy+J2AkSZI0rwzsjvps5x11SZIkDdp0d9Rn35MgkiRJ0jxgUJckSZJayKAuSZIktZBBXZIkSWohg7okSZLUQgZ1SZIkqYUM6pIkSVIL+R71XUiyA/jbIRz6UOCpIRxX+5bXeX7wOs99XuP5wes8PwzrOv9EVS2eaodBvWWSjO7qpfeaO7zO84PXee7zGs8PXuf5oY3X2akvkiRJUgsZ1CVJkqQWMqi3z7XDHoD2Ca/z/OB1nvu8xvOD13l+aN11do66JEmS1ELeUZckSZJayKDeIklWJHkoyViStcMej/qX5Mgkn0/yYJKtSX6rKT8kyW1JHm7WB3e1eW9zrR9KcmZX+euT3NfsuzJJhnFOmlqSBUnuTvKZZttrPMckWZTkU0n+pvk3/TNe57knyb9tfl7fn+TPkxzgdZ79klyX5Mkk93eVzdh1TfLSJH/RlN+V5OhBno9BvSWSLACuAlYCxwFvT3LccEel3bAT+PdV9dPAqcC7muu3Fri9qpYBtzfbNPvOBY4HVgB/3PxvAOBq4BJgWbOs2Jcnop5+C3iwa9trPPf8J+C/VdVrgX9K53p7neeQJEuA3wRGquoEYAGd6+h1nv0+yYuvwUxe14uB/6+qjgU+BPzewM4Eg3qbnAKMVdUjVfU8sBFYNeQxqU9Vtb2q/rr5/AydX+xL6FzD9U219cDq5vMqYGNV/aCqvgmMAackOQI4qKq+Up0HSDZ0tdGQJVkKnA18rKvYazyHJDkIeBPwcYCqer6qvo3XeS5aCByYZCHwY8A2vM6zXlV9Efhfk4pn8rp29/Up4PRB/hXFoN4eS4DHurbHmzLNMs2fwV4H3AUcXlXboRPmgcOaaru63kuaz5PL1Q4fBn4HeKGrzGs8t7wa2AF8opni9LEkL8PrPKdU1ePAHwCPAtuBp6vqv+N1nqtm8rr+qE1V7QSeBl45qIEb1Ntjqv835it5ZpkkLwf+Enh3VX1nuqpTlNU05RqyJG8BnqyqLf02maLMa9x+C4GTgaur6nXAszR/Jt8Fr/Ms1MxRXgUcA7wKeFmSX5uuyRRlXufZb0+u6z695gb19hgHjuzaXkrnz3CaJZK8hE5I/89V9VdN8RPNn9Bo1k825bu63uPN58nlGr5/BvxSkm/RmZr2C0n+DK/xXDMOjFfVXc32p+gEd6/z3PK/A9+sqh1V9ffAXwE/i9d5rprJ6/qjNs20qVfw4qk2M8ag3h6bgWVJjkmyP52HGzYNeUzqUzM/7ePAg1X1h127NgEXNp8vBG7qKj+3eXr8GDoPqny1+ZPcM0lObfq8oKuNhqiq3ltVS6vqaDr/Pv9HVf0aXuM5par+DngsyWuaotOBB/A6zzWPAqcm+bHm+pxO59kir/PcNJPXtbuvX6Xzu2Bgd9QXDqpj7Z6q2plkDXArnafPr6uqrUMelvr3z4B3APcluacpuwxYB9yQ5GI6vxjOAaiqrUluoBMAdgLvqqofNu3eSeep9QOBW5pF7eU1nnv+DfCfm5smjwAX0bmx5XWeI6rqriSfAv6aznW7m863Ur4cr/OsluTPgTcDhyYZB97PzP6c/jjwp0nG6NxJP3eg5+M3k0qSJEnt49QXSZIkqYUM6pIkSVILGdQlSZKkFjKoS5IkSS1kUJckSZJayKAuSdojSVYnOW6G+vpkkl+dib4kaa4wqEuS9tRqYEaCuiTpxQzqkjQPJTk6yd8kWZ/k3iSfar6l8fQkdye5L8l1SV7a1F+X5IGm7h8k+Vngl4DfT3JPkp9McmySzyX5WpK/bsqS5PeT3N/0+S+a/pLkI02fnwUO6xrb65N8IcmWJLdOfPW3JM03fuGRJM1DSY4GvgmcVlVfTnIdnW/h/HXg9Kr6epINdL65cQPwFeC1VVVJFlXVt5N8EvhMVX2q6fMuYF1V3ZjkADo3g1YClwIrgEOBzcBy4GfofPPfCuBwOt8M+K/pfE33F4BVVbWjCfZnVtW/Gvx/FUlqF++oS9L89VhVfbn5/GfA6cA3q+rrTdl64E3Ad4DvAx9L8svAc5M7SvLjwJKquhGgqr5fVc8BpwF/XlU/rKon6ITwNzT9TpRvA/5H09VrgBOA25LcA7wPWDrD5y1Js8LCYQ9AkjQ0ff1Jtap2JjmFTpA/F1gD/MKkatlF812V7+r4AbZW1c/0MzZJmsu8oy5J89dRSSYC8duBzwFHJzm2KXsH8IUkLwdeUVU3A+8GTmr2PwP8OEBVfQcYT7IaIMlLk/wY8EXgXyRZkGQxnTvpX23Kz23KjwB+vunzIWDxxLiSvCTJ8QM5e0lqOYO6JM1fDwIXJrkXOAT4EHAR8F+S3Ae8AFxDJ4x/pqn3BeDfNu03Av9H8/DpT9IJ9r/Z1PufwD8BbgTuBb5GZ3rL71TV3zXlDwP3AVc3/VJVzwO/Cvxekq8B9wA/O8j/CJLUVj5MKknzUPMw6Weq6oRhj0WSNDXvqEuSJEkt5B11SZIkqYW8oy5JkiS1kEFdkiRJaiGDuiRJktRCBnVJkiSphQzqkiRJUgsZ1CVJkqQW+v8BhoQ9ImnfWa0AAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 864x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, axes = plt.subplots(1, figsize=(12, 5))\n",
    "plt.plot(samplespace, p, '.', )\n",
    "axes.set_xlabel('postcode')\n",
    "axes.set_ylabel('probability')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Exercise incorporating more prior knowledge\n",
    "\n",
    "Now try to incorporate the additional prior knowledge that 40% of all mail\n",
    "goes to the following CBD postcodes:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "attributes": {
     "classes": [],
     "id": "",
     "n": "40"
    }
   },
   "outputs": [],
   "source": [
    "CBD_POSTCODES = {2000, 2001, 3000, 3001, 4000, 4001, 5000, 5001, 6000, 6001}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Solution hint:**\n",
    "1. Define a new feature function `in_cbd(postcode)` and append this to your list of features.\n",
    "2. Add an additional value (0.4) to your array of constraint values of feature expectations.\n",
    "3. Re-create your model passing in your new features.\n",
    "4. Re-fit your model passing in the new constraints."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*For solutions, see `solutions/representing_prior_information.py`*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "ename": "NameError",
     "evalue": "name 'k' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "File \u001b[1;32m~\\Documents\\Prework\\Python\\QuantFinance\\notebooks\\solutions\\representing_prior_information.py:7\u001b[0m\n\u001b[0;32m      3\u001b[0m     \u001b[38;5;28;01mreturn\u001b[39;00m [postcode \u001b[38;5;129;01min\u001b[39;00m CBD_POSTCODES \u001b[38;5;28;01mfor\u001b[39;00m postcode \u001b[38;5;129;01min\u001b[39;00m postcodes]\n\u001b[0;32m      5\u001b[0m features2 \u001b[38;5;241m=\u001b[39m features \u001b[38;5;241m+\u001b[39m [in_cbd]\n\u001b[1;32m----> 7\u001b[0m k2 \u001b[38;5;241m=\u001b[39m np\u001b[38;5;241m.\u001b[39mc_[\u001b[43mk\u001b[49m, \u001b[38;5;241m0.25\u001b[39m]\n\u001b[0;32m      8\u001b[0m k2\n\u001b[0;32m     10\u001b[0m \u001b[38;5;28;01massert\u001b[39;00m \u001b[38;5;28mlen\u001b[39m(features2) \u001b[38;5;241m==\u001b[39m k2\u001b[38;5;241m.\u001b[39mshape[\u001b[38;5;241m1\u001b[39m]\n",
      "\u001b[1;31mNameError\u001b[0m: name 'k' is not defined"
     ]
    }
   ],
   "source": [
    "%run -i solutions/representing_prior_information.py"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### More prior knowledge: CBD, inner suburbs, outer suburbs, regional centres\n",
    "\n",
    "Here is an example of how to incorporate this extra information:\n",
    "\n",
    "- Within each state xxxx, 80% of mail goes to x0xx and x1xx (metropolitan city areas and suburbs)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "attributes": {
     "classes": [],
     "id": "",
     "n": "56"
    }
   },
   "outputs": [],
   "source": [
    "def which_ring(postcodes):\n",
    "    \"\"\"\n",
    "    Returns\n",
    "    -------\n",
    "    0 if postcode is x0xx\n",
    "    100 if postcode is x1xx\n",
    "    200 if postcode is x2xx\n",
    "    ... otherwise\n",
    "    \"\"\"\n",
    "    return [postcode % 1000 - postcode % 100 for postcode in postcodes]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "attributes": {
     "classes": [],
     "id": "",
     "n": "57"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[200, 800, 900, 0, 0, 0, 100]"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "which_ring([1234, 800, 2900, 3000, 2001, 2099, 3122])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "attributes": {
     "classes": [],
     "id": "",
     "n": "58"
    }
   },
   "outputs": [],
   "source": [
    "def in_city_metropolitan_area(postcodes):\n",
    "    return [ring == 0 or ring == 100 for ring in which_ring(postcodes)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {
    "attributes": {
     "classes": [],
     "id": "",
     "n": "59"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[True, True, True]"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "in_city_metropolitan_area([3136, 3122, 2001])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {
    "attributes": {
     "classes": [],
     "id": "",
     "n": "60"
    }
   },
   "outputs": [],
   "source": [
    "features3 = features2 + [in_city_metropolitan_area]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {
    "attributes": {
     "classes": [],
     "id": "",
     "n": "61"
    }
   },
   "outputs": [],
   "source": [
    "k3 = np.c_[k2, 0.8]  # k2 comes from your solution to the previous exercises"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[1.   , 0.32 , 0.252, 0.201, 0.071, 0.108, 0.021, 0.01 , 0.016,\n",
       "        0.25 , 0.8  ]])"
      ]
     },
     "execution_count": 40,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "k3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {
    "attributes": {
     "classes": [],
     "id": "",
     "n": "62"
    }
   },
   "outputs": [
    {
     "ename": "NameError",
     "evalue": "name 'MinDivergenceModel' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "Cell \u001b[1;32mIn[29], line 1\u001b[0m\n\u001b[1;32m----> 1\u001b[0m model3 \u001b[38;5;241m=\u001b[39m \u001b[43mMinDivergenceModel\u001b[49m(features3, samplespace, vectorized\u001b[38;5;241m=\u001b[39m\u001b[38;5;28;01mTrue\u001b[39;00m)\n",
      "\u001b[1;31mNameError\u001b[0m: name 'MinDivergenceModel' is not defined"
     ]
    }
   ],
   "source": [
    "model3 = MinDivergenceModel(features3, samplespace, vectorized=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {
    "attributes": {
     "classes": [],
     "id": "",
     "n": "63"
    }
   },
   "outputs": [],
   "source": [
    "model3.fit(k3);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {
    "attributes": {
     "classes": [],
     "id": "",
     "n": "64"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 43,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.allclose(model3.expectations(), k3, atol=1e-5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {
    "attributes": {
     "classes": [],
     "id": "",
     "n": "65"
    }
   },
   "outputs": [],
   "source": [
    "p3 = model3.probdist()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {
    "attributes": {
     "classes": [],
     "id": "",
     "n": "66"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, 'probability')"
      ]
     },
     "execution_count": 45,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAtkAAAE9CAYAAADecH6XAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAltUlEQVR4nO3df5icdXnv8fdN0ghC3YDEqgkx4EYEY4tlTazlUCzlEH6s0GrbQNW2UFI9pWp/R/TIwVKTXlVrLRxojlBAK5SiUgIBSq0aqTQkkSjBCAkRYQFNqBLFYmOS+/wxs7gu+2My+8w+zzz7fl3XXLvz3Xm+c0+ehPnw3Xu+T2QmkiRJkoqzX9kFSJIkSXVjyJYkSZIKZsiWJEmSCmbIliRJkgpmyJYkSZIKZsiWJEmSCja97AI64dBDD8158+aVXYYkSZJqbsOGDU9k5qzh47UK2RHRD/T39vayfv36ssuRJElSzUXEN0Yar1W7SGauysylPT09ZZciSZKkKaxWITsi+iNi5c6dO8suRZIkSVNYrUK2K9mSJEmqglqFbFeyJUmSVAW1CtmuZEuSJKkKahWyXcmWJElSFdQqZLuSLUmSpCqoVciWJEmSqqBWIdt2EXXaJ9Y+zJuvWMsn1j5cdimSJKnCanXFx8xcBazq6+s7r+xaVD+fWPswF3z6XgC+sOUJAM5eNLfMkiRJUkXVaiVb6qRbNz0+5n1JkqRBtQrZtouok05Z8KIx70uSJA2yXURq0WBryK2bHueUBS+aMq0in1j78JR7zZIkTVStQrbUaWcvmjulgqZ96JIktadW7SKSimUfuiRJ7alVyLYnWyrWVO1Dd6tGSdJE1apdxJ5sqVhTsQ/dFhlJUhFqFbIlFW+q9aGP1CIzlV6/JKkYtWoXkaSJmqotMpKkYrmSLUlDTMUWGUlS8QzZkjTMVGuRkSQVr1btIu4uIkmSpCqoVcjOzFWZubSnp6fsUiRJkjSF1SpkS5IkSVVgyJYkTUledEhSJ/nBR0nSlONFhyR1mivZkqQpZ6SLDklSkboiZEfEmRHx/yLinyPif5ZdjySpu3nRIUmd1vF2kYi4Ejgd2J6ZC4aMLwb+BpgGfDQzV4w2R2beCNwYEQcDHwD+paNFS5JqzYsOSeq0yMzOPkHE8cBTwDWDITsipgEPACcBA8A64CwagXv5sCnOycztzeM+CPxDZn5prOfs6+vL9evXF/o6JEmSpOEiYkNm9g0f7/hKdmauiYh5w4YXAlszc1uzuOuAMzJzOY1V7x8TEQGsAG4dL2BLkvbdJ9Y+7KquJBWorN1FZgOPDLk/ACwa4/G/D/wS0BMRvZl5+fAHRMRSYCnA3Lm+QUhSq9xpQ5KKV9YHH2OEsVH7VjLzI5l5bGa+daSA3XzMyszsy8y+WbNmFVaoJNWdO21IUvHKCtkDwGFD7s8BHpvopBHRHxErd+7cOdGpJGnKcKcNSSpeWe0i64D5EXE48CiwBDh7opNm5ipgVV9f33kTnUuSpgp32pCk4k3GFn7XAicAh0bEAHBhZl4REecDt9PYUeTKzLyvgOfqB/p7e3snOpUkTSlnL5pruJakAnV8C78yuIWfJEmSJsNoW/h1xRUfW2VPtiRJkqqgViE7M1dl5tKenp6yS5EkSdIUVquQ7Uq2JEmSqqBWIduVbEmSJFVBrUK2JEmSVAW1Ctm2i0iSJKkKahWybReRJElSFdQqZEuSJElVUKuQbbuIJEmSqqBWIdt2EUmSJFVBrUK2JEmSVAWGbEmSJKlgtQrZ9mRLkiSpCmoVsu3JliRJUhXUKmRLkiRJVWDIliRJkgpmyJYkSZIKZsiWJEmSClarkO3uIpIkSaqCWoVsdxeRJElSFdQqZEuSJElVYMiWJEmSCmbIliRJkgpmyJYkSZIKVvmQHRFHRcTlEXFDRLyt7HokSZKk8XQ0ZEfElRGxPSI2DRtfHBH3R8TWiFg21hyZuTkz3wr8GtDXyXolSZKkInR6JfsqYPHQgYiYBlwKnAIcDZwVEUdHxCsj4uZhtxc0j3k9cCfwmQ7XK0mSJE3Y9E5OnplrImLesOGFwNbM3AYQEdcBZ2TmcuD0Uea5CbgpIm4BPtHBkiVJkqQJ62jIHsVs4JEh9weARaM9OCJOAH4FeA6weozHLQWWAsydO7eAMqVne8sVa1mz5Yl9Pu6A6fux+eJTOlCRJEmqojI++BgjjOVoD87Mz2Xm2zPzdzPz0jEetzIz+zKzb9asWYUUKg3VbsAGeHr3Xo56z60FVyRJkqqqjJA9ABw25P4c4LEiJo6I/ohYuXPnziKmk37M3Q99e0LHP717b0GVSJKkqiujXWQdMD8iDgceBZYAZ5dQh7RPFs47pO2VbGi0jHSjkz74Obbs+P4+HzfzgOlsvPDkDlQkSVL1dXoLv2uBu4AjI2IgIs7NzN3A+cDtwGbg+sy8r4jny8xVmbm0p6eniOmkH3PNuYs4fv6hbR3brT3Z7QZsgCef3s0xF91ecEWSJHWHTu8uctYo46sZ40OM7YqIfqC/t7e36KkloBG0p5IHn2gvYA968undBVUiSVJ36c7fX4/ClWypWC899MAJHT/zgDI60iRJKp/vgJJGdccfnTAle7JfffEd7Hhq1z4fN2fm/ty57MQOVCRJ6jaROerueV1nSLvIeVu2bCm7HEldqN2APcigLUlTS0RsyMy+4eO2i0jSEBMJ2ACPPvmDgiqRJHWzWoVs98mWNFGzDpoxoeNnz9y/oEokSd2sViHblWxJE7XuPSe1HbRtFZEkDfKDj1KXe+m7bmFP86MVx88/tO1tBo+56PZnttzr5g8tFmHde04quwRJUper1Uq2NNUMDdgAa7Y8wVuuWLvP8wwN2OCFZCRJmqhahWx7sjXV7Blhc6C7H/r2Ps8z0kVjvJCMJEntq1W7SGauAlb19fWdV3Yt0mSYFs8O2gvnHbLP88w8YPqzQrUXklHdvezdq9k10v+pjuPMY17Mh5e8qgMVSaqTWq1kS1PNg8tPY1r86H67PdkbLzz5x0L1VO/JVv21G7ABbtz4GO+87p6CK5JUNy5VSV3uweWnFTKPoVpTSbsBe9DnHthRUCWS6qpWK9n2ZEuSWjFj6K+A2nDCy2YVVImkuqpVyHafbElSKx74i1PbDtr2ZEtqhe0ikqQp6YG/OLXsEiTVmCFbksSC997GU7v27PNxx8zp4cbzj+tARZLU3WrVLiJJ2nftBmyAjQM7OfOSOwuuSJK6nyFbkqa4dgP2oE2PfbegSiSpPmoVst1dRJL23UEzpk3o+AUvfl5BlUhSfdQqZLu7iCTtu03vW9x20LYnW5JG5gcfJUlset/iskuQpFqp1Uq2JEmSVAWGbEmSJKlghmxJkiSpYF0RsiPiwIjYEBGnl12LJEmSNJ6OhuyIuDIitkfEpmHjiyPi/ojYGhHLWpjqz4DrO1OlJEmSVKxO7y5yFXAJcM3gQERMAy4FTgIGgHURcRMwDVg+7PhzgJ8Gvgrs3+Fap7zjVnyGgSd/MOrPZx4wnY0XnjyJFUmSJHWnjobszFwTEfOGDS8EtmbmNoCIuA44IzOXA89qB4mI1wEHAkcDT0fE6szc28m6p6LxAjbAk0/v5piLbjdoS5IkjaOMfbJnA48MuT8ALBrtwZn5boCI+C3gidECdkQsBZYCzJ07t6hap4xHxwnYg558eneHK5EkdcKGb3yHN1z2xbaOff8vv5KzF/neKu2LMkJ2jDCW4x2UmVeN8/OVEfE40D9jxoxj26xtypo9c/9xV7Kh0TIiSeouEwnYABd8+l4Ag7a0D8pITAPAYUPuzwEeK2LizFwFrOrr6zuviPmmkjuXnWhPdge887p7uHFjIX+9Szd/1oHc8UcnlF2GpDb8x7b/nPAct2563JAt7YMyQvY6YH5EHA48CiwBzi5i4ojoB/p7e3uLmG7KuXPZiWWXUCt1CtgAW3Z8n5M++DmDttSFXnPE8yc8xykLXlRAJdLU0ekt/K4F7gKOjIiBiDg3M3cD5wO3A5uB6zPzviKeLzNXZebSnp6eIqaTJuRzD+wou4TCPfjE98suQVIbjn3JwXzyba9t+3h7sqV91+ndRc4aZXw1sLro53MlW1Vywstm1WolG+Clhx5YdgmS2nTsSw7moRWnlV2GNGXU6lNs9mSrSj685FUAtQnarfZkn3nJnWwc2Nn5gibB8fMP5ZpzR938SJKkUUXmuBt7dI0hK9nnbdmypexypCmnTgF7kEFbkjSWiNiQmX3Dxzvakz3Z7MmWyrXpse+WXULh7n7o22WXIEnqQrVqF5FUrgUvfl7tVrIXzjtk3Me0csXUbvHW449g2alHlV2GJHW9Wq1kR0R/RKzcubNeb/JSt7jx/OM4Zk59fpPUSqtInQI2wOVrtrFi9eayy5CkrtdST3ZEfBK4Erh1tMuaV0lfX1+uX7++7DIkTQGHL7tl/EvWdpl5z38un/uT15VdhiR1hYn2ZF9G44IxWyJiRUS8vNDqJKlLzZ65f9klFG7xK15YdgmS1PVaCtmZ+a+Z+RvAzwIPAXdExBcj4rcj4ic6WeC+sF1E0mS7c9mJzKlR0LYnW5KK0fIWfhHxfOBNwJuBx4B/AI4DXpmZJ3SqwHbYLiJJkqTJMFq7SEu7i0TEp4CXAx8D+jPz8eaP/jEiTLOSJEnSEK1u4ffR5qXQnxERz8nM/x4puUuSJElTWash+2Jg9bCxu2j0aEuSVHtHvedWnt5dvQ22DpoxjU3vW1x2GZKGGfODjxHxwog4FjggIl4VET/bvJ0APHcyCtwXfvBRktQJVQ3YAE/t2sOC995WdhmShhlvJftk4LeAOcCHhox/D7igQzW1LTNXAav6+vrOK7sWSVJ9VDVgD3pq156yS5A0zJghOzOvBq6OiDdk5icnqSZJkirlgOn7VTpoHzRjWtklSBpmzJAdEW/KzI8D8yLiD4f/PDM/NMJhkiTVyuaLT6lsy4g92VI1jdcucmDz60GdLkSSpCrbfPEpZZcgqYuM1y7yd82vF01OOZKkbnDMRbfz5NO7yy7jWaYFPLj8tLLLkKRx20U+MtbPM/PtxZYzMRHRD/T39vaWXYok1VZVAzbAnoSXvusWg7ak0o3XLrJhUqooiLuLSFLnVTVgD9qTZVcgSa3tLiJJ0jNmHjC90kF7WpRdgSSNfzGaDze/roqIm4bfJqVCSVKlbLzwZGYe0OoFgyeXPdmSqmK8/0p+rPn1A50uRJLUPTZeeHLZJUhSpY3XLrKh+fXzETEDeDmQwP2ZuWsS6pMkSZK6zpjtIoMi4jTgQeAjwCXA1oiYlA1DI+KEiPhCRFweESdMxnNKkiRJE9FSyAY+CLwuM0/IzF8AXgf89XgHRcSVEbE9IjYNG18cEfdHxNaIWDbONAk8BewPDLRYryRJklSaVj+5sj0ztw65vw3Y3sJxV9FY+b5mcCAipgGXAifRCM3rmh+inAYsH3b8OcAXmu0qPwV8CPiNFmuWJEmSSjHexWh+pfntfRGxGriexsryrwLrxps8M9dExLxhwwuBrZm5rfkc1wFnZOZy4PQxpvsO8JzxnlOSJEkq23gr2f1Dvv8W8AvN73cAB7f5nLOBR4bcHwAWjfbgZtA/GZhJY1V8tMctBZYCzJ07t83SJEmSpIkbb3eR3+7Ac450mYBRr8+VmZ8CPjXepJm5MiIeB/pnzJhx7ATqkyRJkiakpZ7siNgfOBd4BY0PIAKQmee08ZwDwGFD7s8BHmtjnmfxsuqSJEmqglY/+Pgx4Gs02jbeR+PDh5vbfM51wPyIOBx4FFgCnN3mXD8mIvqB/t7e3iKmkyRpSuu94BZ27y27imc7aMY0Nr1vcdllSGNqdQu/3sz838D3M/Nq4DTgleMdFBHXAncBR0bEQEScm5m7gfOB22kE9esz8772yv9xmbkqM5f29PQUMZ0kSVNWVQM2wFO79rDgvbeVXYY0plZXsn/Y/PpkRCwAvgnMG++gzDxrlPHVwOoWn7tlrmRLklSMqgbsQU/t2lN2CdKYWl3JXhkRBwP/G7gJ+Crwlx2rqk2uZEuSVIzprSaEkhw0Y1rZJUhjaumfUGZ+NDO/k5mfz8wjMvMFmfl3nS5uX0VEf0Ss3LlzZ9mlSJLU1ba+/7TKBm17stUNInPU3fN+9KCI5wP/B/h5GtvtfQH488z8z45W16a+vr5cv3592WVIkiSp5iJiQ2b2DR9v9f9Rr6NxGfU3AG8EngD+sbjyJEmSpPpoNWQfkpl/nplfb94upnEFxkqxXUSSJElV0GrI/mxELImI/Zq3XwNu6WRh7fCDj5IkSaqCMbfwi4jv0ejBDuAPgY83f7Qf8BRwYUerkyRJkrrQmCE7M39ysgopgvtkS5IkqQpa3pwnIl4fER9o3k7vZFHtsl1EkiRJVdBSyI6IFcA7aFyE5qvAO5pjkiRJkoZp9bLqpwLHZOZegIi4GrgHWNapwiRJkqRutS/Xcpo55PtK9mO4hZ8kSZKqoNWQ/X7gnoi4qrmKvaE5Vin2ZEuSJKkKxm0XiYj9gL3Aa4BX09jO788y85sdrk2SJEnqSuOG7MzcGxHnZ+b1wE2TUJMkSZLU1VptF7kjIv44Ig6LiEMGbx2tTJIkSepSre4ucg6NKz/+r2HjRxRbjiRJktT9Wl3JPhq4FPgysBH4W+AVHaqpbe4uIkmSpCpoNWRfDRwFfIRGwD6qOVYp7i4iSZKkKmi1XeTIzPyZIfc/GxFf7kRBkiRJUrdrdSX7noh4zeCdiFgE/HtnSpIkSZK6W6sr2YuAt0TEw837c4HNEXEvkJn50x2pTpIkSepCrYbsxR2tQpIkSaqRlkJ2Zn6j04VIkiRJddHqSnZpmpd1/3PgecD6zKzcriZT1VuuWMuaLU/s83HT94Ot7z+tAxVJkvQjK1Zv5vI12/b5uAC+vqI736fafW+eFvDg8u58zVXV6gcf2xIRV0bE9ojYNGx8cUTcHxFbI2LZONOcAcwGfggMdKpW7Zt2/xED7N4LvRfcUnBFkiT9SLsBGxpX3zt8Wfe9T03kvXlPwkvf1X2vuco6vZJ9FXAJcM3gQERMo3Fhm5NohOZ1EXETMA1YPuz4c4Ajgbsy8+8i4gbgMx2uWS24+6FvT+j43XsLKkSSpBHcdt83J3R8FlTHZJroe/OebnzRFdbRlezMXAMMP+MLga2ZuS0zdwHXAWdk5r2Zefqw23YaQfw7zWP3jPZcEbE0ItZHxPodO3Z04uVoiIXzDpnQ8dM7+jdPkjTVLX7FCyd0fBRUx2Sa6HvztG580RVWRtSZDTwy5P5Ac2w0nwJOjoi/BdaM9qDMXJmZfZnZN2vWrGIq1aiuOXcRx88/tK1j7cmWJHXaslOP4q3HH9HWsd3akz2R92Z7sosXmZ393UBEzANuzswFzfu/Cpycmb/TvP9mYGFm/n4Bz9UP9Pf29p63ZcuWiU4nSZIkjSkiNmRm3/DxMlayB4DDhtyfAzxWxMSZuSozl/b09BQxnSRJktSWMkL2OmB+RBweETOAJcBNRUwcEf0RsXLnzp1FTCdJkiS1pdNb+F0L3AUcGREDEXFuZu4GzgduBzYD12fmfUU8nyvZkiRJqoKObuGXmWeNMr4aWF308w3pyS56akmSJKlltdpIzZVsSZIkVUGtQrY92ZIkSaqCWoVsV7IlSZJUBbUK2ZIkSVIV1Cpk2y4iSZKkKqhVyLZdRJIkSVVQq5AtSZIkVUGtQrbtIpIkSaqCWoVs20UkSZJUBbUK2ZIkSVIVGLIlSZKkgtUqZNuTLUmSpCqoVci2J1uSJElVUKuQLUmSJFWBIVuSJEkqmCFbkiRJKpghW5IkSSpYrUK2u4tIkiSpCmoVst1dRJIkSVVQq5AtSZIkVYEhW5IkSSqYIVuSJEkqmCFbkiRJKtj0sgsYT0T8D+A3aNR6dGa+tuSSJEmSpDF1dCU7Iq6MiO0RsWnY+OKIuD8itkbEsrHmyMwvZOZbgZuBqztZryRJklSETq9kXwVcAlwzOBAR04BLgZOAAWBdRNwETAOWDzv+nMzc3vz+bOB3OlyvJEmSNGEdDdmZuSYi5g0bXghszcxtABFxHXBGZi4HTh9pnoiYC+zMzO92st7JdviyW8iyixjBrINmsO49Jz1z/8xL7mTjwOgX+DloxjQ2vW/xuPN+Yu3DXPDpewupsWzD/4wklWvDN77DGy77YlvHfvJtr+XYlxxccEWd95Yr1rJmyxP7fNx+wLYVpxVf0CRo9zVP3w+2vr+113zEslvYu8/PUD0HTN+PzRefMu7j6vTePPOA6Wy88OSyy3hGGR98nA08MuT+QHNsLOcCfz/WAyJiaUSsj4j1O3bsmGCJnVfVgA2w46ldvPriO4DxAzbAU7v2sOC9t435mDr9I4Yf/zOSVK6JBGyAN1z2RTZ84zsFVtR57YZNgL00gmS3mchr3r0Xei8Y/zXXJWADPL17L0e959YxH1O39+Ynn97NMRfdXnYZzygjZMcIY2Pmzcy8MDPH/C9oZq7MzL7M7Js1a9aECpwMVQ3Yg3Y8tQuATY+19suDp3btGfPnt256fMI1Vc3gn5Gkcv3Htv+sxByT6e6Hvj2h47sxSE70Ne9u4UV345/LWJ4e50XX8b35yad3l13CM8oI2QPAYUPuzwEeK2LiiOiPiJU7d4698loFI/2fRpXMOmgGAAte/LyWHn/QjGlj/vyUBS+acE1VM/hnJKlcrzni+ZWYYzItnHfIhI7vxv17J/qap7fworvxz2UsB4zzouv43jzzgOpsnFfG36d1wPyIODwiZgBLgJtKqKNUX19xWmWD9tB+4xvPP45j5vSM+fhWerLPXjSX9//yKwursWz2ZEvVcexLDuaTb2t/d9du7Mm+5txFHD//0LaO7dae7Im85lZ7sretOK02QbuVnuy6vTdXrSc7MjvXuBAR1wInAIcC3wIuzMwrIuJU4MM0dhS5MjP/osjn7evry/Xr1xc5pSRJkvQsEbEhM/uGj3d6d5GzRhlfDawu+vkioh/o7+3tLXpqSZIkqWV1+a0IAJm5KjOX9vSM3d4gSZIkdVKtQrYkSZJUBbUK2d20u4gkSZLqq1Yh23YRSZIkVUGtQrYr2ZIkSaqCWoVsV7IlSZJUBbUK2ZIkSVIVGLIlSZKkgtUqZNuTLUmSpCqoVci2J1uSJElVUKuQLUmSJFWBIVuSJEkqWK1Ctj3ZkiRJqoJahWx7siVJklQFtQrZkiRJUhUYsiVJkqSCGbIlSZKkghmyJUmSpILVKmS7u4gkSZKqoFYh291FJEmSVAW1CtmSJElSFRiyJUmSpIIZsiVJkqSCGbIlSZKkgk0vu4DxRMRc4BLgCeCBzFxRckmSJEnSmDq6kh0RV0bE9ojYNGx8cUTcHxFbI2LZONO8DLglM88Bju5YsZIkSVJBOt0uchWweOhAREwDLgVOoRGaz4qIoyPilRFx87DbC4B7gCUR8W/AZztcryRJkjRhHW0Xycw1ETFv2PBCYGtmbgOIiOuAMzJzOXD68Dki4o+BC5tz3QD8fSdrliRJkiaqjA8+zgYeGXJ/oDk2mtuAt0fE5cBDoz0oIpZGxPqIWL9jx45CCpUkSZLaUcYHH2OEsRztwZm5CXjjeJNm5kpgJUBfX9+o80mSJEmdVsZK9gBw2JD7c4DHipg4IvojYuXOnTuLmE6SJElqSxkhex0wPyIOj4gZwBLgphLqkCRJkjqi01v4XQvcBRwZEQMRcW5m7gbOB24HNgPXZ+Z9RTxfZq7KzKU9PT1FTCdJkiS1pdO7i5w1yvhqYHXRzxcR/UB/b29v0VNLkiRJLavVZdVdyZYkSVIV1Cpk+8FHSZIkVUGtQrYr2ZIkSaqCWoVsSZIkqQpqFbJtF5EkSVIV1Cpk2y4iSZKkKqhVyJYkSZKqoFYh23YRSZIkVUGtQrbtIpIkSaqCWoVsSZIkqQoM2ZIkSVLBDNmSJElSwWoVsv3goyRJkqqgViHbDz5KkiSpCmoVsiVJkqQqmF52AXWxYvVmLl+zbZ+P2w/YtuK04guSJElSaVzJLkC7ARtgL3DEsluKLUiSJEmlMmQX4Lb7vjmh4/cWVIckSZKqoVYhu6zdRRa/4oUTOr5WJ0GSJEn1yndl7S6y7NSjeOvxR7R1rD3ZkiRJ9ROZWXYNhevr68v169eXXYYkSZJqLiI2ZGbf8PFarWRLkiRJVWDIliRJkgpmyJYkSZIKVvmQHRFHR8T1EXFZRLyx7HokSZKk8XQ0ZEfElRGxPSI2DRtfHBH3R8TWiFg2zjSnAH+bmW8D3tKxYiVJkqSCdPqy6lcBlwDXDA5ExDTgUuAkYABYFxE3AdOA5cOOPwf4GHBhRLweeH6H65UkSZImrKMhOzPXRMS8YcMLga2ZuQ0gIq4DzsjM5cDpo0z1e81w/qmOFStJkiQVpNMr2SOZDTwy5P4AsGi0BzdD+gXAgcBfjfG4pcBSgLlz5xZRpyRJktSWMkJ2jDA26hVxMvMhmuF5LJm5MiIeB/pnzJhxbPvlSZIkSRNTRsgeAA4bcn8O8FgRE2fmKmBVRPxyRHyjiDn30aHAEyU8ryaX53lq8DxPDZ7n+vMcTw1lnueXjDRYRsheB8yPiMOBR4ElwNlFPkFmzipyvlZFxPqRLqupevE8Tw2e56nB81x/nuOpoYrnudNb+F0L3AUcGREDEXFuZu4GzgduBzYD12fmfZ2sQ5IkSZpMnd5d5KxRxlcDqzv53JIkSVJZKn/Fxy6zsuwCNCk8z1OD53lq8DzXn+d4aqjceY7MUTf2kCRJktQGV7IlSZKkghmyCxIRiyPi/ojYGhHLyq5HrYuIwyLisxGxOSLui4h3NMcPiYg7ImJL8+vBQ455V/Nc3x8RJw8ZPzYi7m3+7CMRMdK+8CpRREyLiHsi4ubmfc9zzUTEzIi4ISK+1vx3/XOe53qJiD9o/vd6U0RcGxH7e467X0RcGRHbI2LTkLHCzmtEPCci/rE5vnaEq5IXypBdgOYl3y8FTgGOBs6KiKPLrUr7YDfwR5l5FPAa4Pea528Z8JnMnA98pnmf5s+WAK8AFgP/t/l3AOAyGhdPmt+8LZ7MF6KWvIPGzkaDPM/18zfAbZn5cuBnaJxvz3NNRMRs4O1AX2YuAKbROIee4+53Fc8+B0We13OB72RmL/DXwF927JVgyC7KQmBrZm7LzF3AdcAZJdekFmXm45n5peb336Pxhjybxjm8uvmwq4Ezm9+fAVyXmf+dmV8HtgILI+JFwPMy865sfNjhmiHHqAIiYg5wGvDRIcOe5xqJiOcBxwNXAGTmrsx8Es9z3UwHDoiI6cBzaVzUznPc5TJzDfDtYcNFntehc90AnNjJ314YsosxG3hkyP2B5pi6TPNXR68C1gI/lZmPQyOIAy9oPmy08z27+f3wcVXHh4E/BfYOGfM818sRwA7g75ttQR+NiAPxPNdGZj4KfAB4GHgc2JmZ/4LnuK6KPK/PHNO8bstO4PmdKtyQXYyR/i/IbVu6TEQcBHwSeGdmfnesh44wlmOMqwIi4nRge2ZuaPWQEcY8z9U3HfhZ4LLMfBXwfZq/Xh6F57nLNHtyzwAOB14MHBgRbxrrkBHGPMfdr53zOqnn3JBdjAHgsCH359D41ZW6RET8BI2A/Q+Z+anm8Leav3ai+XV7c3y08z3Q/H74uKrh54HXR8RDNFq6fjEiPo7nuW4GgIHMXNu8fwON0O15ro9fAr6emTsy84fAp4DX4jmuqyLP6zPHNFuNenh2e0phDNnFWAfMj4jDI2IGjUb8m0quSS1q9mNdAWzOzA8N+dFNwG82v/9N4J+HjC9pfkr5cBofqri7+Wus70XEa5pzvmXIMSpZZr4rM+dk5jwa/0b/LTPfhOe5VjLzm8AjEXFkc+hE4Kt4nuvkYeA1EfHc5rk5kcZnaTzH9VTkeR061xtpvA907rcXmemtgBtwKvAA8CDw7rLr8bZP5+44Gr8u+gqwsXk7lUaf1meALc2vhww55t3Nc30/cMqQ8T5gU/Nnl9C84JO3at2AE4Cbm997nmt2A44B1jf/Td8IHOx5rtcNuAj4WvP8fAx4jue4+2/AtTT67H9IY9X53CLPK7A/8E80PiR5N3BEJ1+PV3yUJEmSCma7iCRJklQwQ7YkSZJUMEO2JEmSVDBDtiRJklQwQ7YkSZJUMEO2JE1BEXFmRBxd0FxXRcQbi5hLkurCkC1JU9OZQCEhW5L0bIZsSeoyETEvIr4WEVdHxFci4obm1e9OjIh7IuLeiLgyIp7TfPyKiPhq87EfiIjXAq8H/ioiNkbESyOiNyL+NSK+HBFfao5FRPxVRGxqzvnrzfkiIi5pznkL8IIhtR0bEZ+PiA0Rcfvg5ZAlaarxYjSS1GUiYh7wdeC4zPz3iLgS2Ab8LnBiZj4QEdcAXwKuAe4CXp6ZGREzM/PJiLiKxlUvb2jOuRZYkZmfjoj9aSzCnAK8FVgMHAqsAxYBPwe8rTn+UzQuW/47NC5d/HngjMzc0QzlJ2fmOZ3/U5GkanElW5K60yOZ+e/N7z8OnAh8PTMfaI5dDRwPfBf4AfDRiPgV4L+GTxQRPwnMzsxPA2TmDzLzv4DjgGszc09mfotGgH51c97B8ceAf2tOdSSwALgjIjYC7wHmFPy6JakrTC+7AElSW1r6NWRm7o6IhTRC+BLgfOAXhz0sRjl8tPHRnj+A+zLz51qpTZLqzJVsSepOcyNiMMyeBfwrMC8ieptjbwY+HxEHAT2ZuRp4J3BM8+ffA34SIDO/CwxExJkAEfGciHgusAb49YiYFhGzaKxg390cX9IcfxHwuuac9wOzBuuKiJ+IiFd05NVLUsUZsiWpO20GfjMivgIcAvw18NvAP0XEvcBe4HIaQfrm5uM+D/xB8/jrgD9pflDypTRC+dubj/si8ELg08BXgC/TaAn508z8ZnN8C3AvcFlzXjJzF/BG4C8j4svARuC1nfxDkKSq8oOPktRlmh98vDkzF5RdiyRpZK5kS5IkSQVzJVuSJEkqmCvZkiRJUsEM2ZIkSVLBDNmSJElSwQzZkiRJUsEM2ZIkSVLBDNmSJElSwf4/3iA4tx/h7kEAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 864x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, axes = plt.subplots(1, figsize=(12, 5))\n",
    "plt.semilogy(samplespace, p3, '.', )\n",
    "axes.set_xlabel('postcode')\n",
    "axes.set_ylabel('probability')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Conclusion\n",
    "This prior reflects **precisely** the information we put into the model.\n",
    "\n",
    "No less:\n",
    "- all constraints we placed on it are satisfied;\n",
    "\n",
    "No more:\n",
    "- no additional information is reflected / assumed which we didn't explicitly add. It is as flat as possible (maximal entropy) subject to our constraints."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {
    "attributes": {
     "classes": [],
     "id": "",
     "n": "67"
    }
   },
   "outputs": [],
   "source": [
    "np.save('postcode_prior3.npy', p3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Aside: Demonstration that we cannot tweak `model` to be equivalent to `model2` by adding one constraint (`in_cbd`) and then minimizing KL divergence from `model1`.\n",
    "\n",
    "Let's try it ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([-22.32221447, -22.32221447, -22.32221447, ..., -15.27472916,\n",
       "       -15.27472916, -15.27472916])"
      ]
     },
     "execution_count": 47,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model.log_probdist()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([   0,    1,    2, ..., 9997, 9998, 9999], dtype=uint16)"
      ]
     },
     "execution_count": 48,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model.samplespace"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {},
   "outputs": [],
   "source": [
    "def prior_log_pdf(x):\n",
    "    return model.log_probdist()[x]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {
    "attributes": {
     "classes": [],
     "id": "",
     "n": "125"
    }
   },
   "outputs": [],
   "source": [
    "model4 = MinDivergenceModel([in_cbd], samplespace, prior_log_pdf=prior_log_pdf, vectorized=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {
    "attributes": {
     "classes": [],
     "id": "",
     "n": "126"
    }
   },
   "outputs": [],
   "source": [
    "k4 = np.array([0.25], ndmin=2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {
    "attributes": {
     "classes": [],
     "id": "",
     "n": "129"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "MinDivergenceModel(features=[<function in_cbd at 0x7fe5b4252ee0>],\n",
       "                   prior_log_pdf=<function prior_log_pdf at 0x7fe5b51bcca0>,\n",
       "                   samplespace=array([   0,    1,    2, ..., 9997, 9998, 9999], dtype=uint16),\n",
       "                   vectorized=True)"
      ]
     },
     "execution_count": 52,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model4.fit(k4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {
    "attributes": {
     "classes": [],
     "id": "",
     "n": "130"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.24999999])"
      ]
     },
     "execution_count": 53,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model4.expectations()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {
    "attributes": {
     "classes": [],
     "id": "",
     "n": "131"
    }
   },
   "outputs": [],
   "source": [
    "p4 = model4.probdist()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {
    "attributes": {
     "classes": [],
     "id": "",
     "n": "132"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, 'probability')"
      ]
     },
     "execution_count": 55,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAtkAAAFCCAYAAAAkBMGDAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAb/klEQVR4nO3dfbRldXkf8O/jEERBRYFEA4wDSFEkjcYrGGstCWU5iBGrNgGTmEayiDY0JrFNR2N0NTWBrMS8anVRJYBWCPElpYJaNAoaDc6Mr0MIgqgwogESUTExCDz9456xd93MvXOZ2WfOPed+Pmuddc/+nb33ec78uDNf9nn23tXdAQAAhvOASRcAAACzRsgGAICBCdkAADAwIRsAAAYmZAMAwMCEbAAAGJiQDQAAAxOyAQBgYKs+ZFfVkVX15qp6+6RrAQCAlRhryK6q86vqtqratmh8Y1VdX1U3VtWm5fbR3Td195njrBMAAIa0z5j3f0GS1yW5aMdAVa1L8vokJyfZnmRzVV2WZF2ScxZt/6Luvm3MNQIAwKDGGrK7++qq2rBo+PgkN3b3TUlSVZckOa27z0nyrCHe9+CDD+4NGxa/LQAADGfr1q13dPchO3tt3Eeyd+bQJLcsWN6e5ISlVq6qg5L8ZpInVtXLR2F8Z+udleSsJFm/fn22bNkyXMUAALBIVX1pqdcmEbJrJ2O91Mrd/XdJXryrnXb3eUnOS5K5ubkl9wcAAOM2iauLbE9y+ILlw5LcOoE6AABgLCYRsjcnObqqjqiqfZOcnuSyCdQBAABjMe5L+F2c5GNJjqmq7VV1Znffk+TsJO9Lcl2SS7v72nHWAQAAe9O4ry5yxhLjVyS5YpzvDQAAk7Lq7/gIAADTRsgGAICBCdkAADAwIRvuh7ddc3N++s3X5G3X3DzpUgCAVWwSN6OBqfS2a27OK9712STJh2+4I0nyghPWT7IkAGCVciQbVug9276y7DIAwA5CNqzQKcc9atllAIAdtIvACu1oDXnPtq/klOMetWZaRd52zc1r7jMDwJ4SsuF+eMEJ69dU0NSHDgC7R7sIsCR96ACwe4RsYElrtQ/dpRoB2FPaRYAlrcU+dC0yAAxByAaWtdb60HfWIrOWPj8Aw9AuArDAWm2RAWBYjmQDLLAWW2QAGJ6QDbDIWmuRWatcAx4YJyEbgDXHCa7AuOnJBmDNcQ14YNyEbADWHCe4AuOmXQSANccJrsC4CdkArMmTAJ3gCoyTkA2wxjkJEGB4erIB1jgnAQIMT8gGWOOcBAgwPO0iAGuckwABhidkA+AkQICBaRcBAICBCdkAADAwIRsAAAYmZAMAwMCEbAAAGJiQDQAAAxOyAQBgYEI2AAAMTMgGAICBCdkAADAwIRsAAAYmZAMAwMCEbAAAGJiQDQAAAxOyAQBgYEI2AAAMTMgGAICBCdkAADAwIRsAAAYmZAMAwMCEbAAAGJiQDQAAAxOyAQBgYKs+ZFfV46rqjVX19qp6yaTrAQCAXRlryK6q86vqtqratmh8Y1VdX1U3VtWm5fbR3dd194uT/HiSuXHWCwAAQxj3kewLkmxcOFBV65K8PskpSY5NckZVHVtVP1BV7170+N7RNs9O8pEkHxhzvQAAsMf2GefOu/vqqtqwaPj4JDd2901JUlWXJDmtu89J8qwl9nNZksuq6vIkbxtjyQAAsMfGGrKXcGiSWxYsb09ywlIrV9WJSZ6b5IFJrlhmvbOSnJUk69evH6BMAADYPZMI2bWTsV5q5e7+UJIP7Wqn3X1ekvOSZG5ubsn9AQDAuE3i6iLbkxy+YPmwJLdOoA4AABiLSYTszUmOrqojqmrfJKcnuWwCdQAAwFiM+xJ+Fyf5WJJjqmp7VZ3Z3fckOTvJ+5Jcl+TS7r52nHUAAMDeNO6ri5yxxPgVWeYkRgAAmGar/o6PAAAwbYRsAAAYmJANAAADE7IBAGBgQjYAAAxMyAYAgIEJ2QAAMLCxXicbGL+jXn557u35508/+uBcdOYJu7Wf41713tx1971JkkMO2DebX3nyUCUCwJrjSDZMsYUBO0muvuGOvPDN19zv/SwM2Ely+11358mvuXKIEgFgTaru3vVaU2Zubq63bNky6TKYQS988zW5+oY7Jl3GRBx9yP658mUnTroMAFg1qmprd8/t7DVHsmGF1nLATpIbbv9WTn7thyZdBgBMBSEbVujjX/z7SZcwcZ+/41uTLgEApoKQDSt0/IZHTLqEiTvq4P0nXQIATAU92XA/rPWWkbVmT67WAsDsW64nW8gGlnTyaz+UG25f2y0igjYAS3HiI7Bb9GDrxQdg97gZDbCkow7ef80fyf72d+7Lhk2XT7qMvebFTz8ym575uEmXATD1HMkGlnTly07M0Yc42XEteePVN+XcK66bdBkAU8+RbGBZa+0GNGvpqPVS3nvtVx3NBthDjmQDLHDIAftOuoSJ2/j4R066BICpJ2QDLLD5lSev6aCtJxtgGNpFABbZ/MqTJ10Ce8G/+LUrcve9038Z2wP2XZdtv7Fx0mUAiziSDcCaMysBO0nuuvveHPeq9066DGARIRuANWdWAvYOd91976RLABYRsgFYc/ZdV5MuYVAH7Ltu0iUAi+jJ5ruedu4Hsv3Oby/5+oEP2iefevUz9mJFAOPxud985sy0jOjJhtWpuqf/L5jF5ubmesuWLZMuY6rsKmDvIGjDbDruVe+diZaDdZV8/pxTJ10GsEZU1dbuntvZa9pFSJJ8eQUBO0nu/Md7xlwJsLfNSsBOkns7OerlbigETJ6QTZLk0AP3W9F6Bz5IhxHMmlkJ2DvMQAcIMAOEbJIkH9l0Ug7bRdDWKgKzadZOmpuxcxqBKaUnG4CZaRnRkw3sTcv1ZPvuHwBXpwAYmHYRAAAYmJANAAADE7IBAGBgQjYAAAxMyAYAgIEJ2QAAMDAhGwAABiZkAwDAwIRsAAAYmJANAAADc1t1AFgDtn7pa3neGz466TIGccC+67LtNzZOugxYliPZADDjZilgJ8ldd9+b41713kmXAcsSsgFgxv3VTX836RIGd9fd9066BFjWikJ2Vb2jqk6tKqEcAKbMU448aNIlDO6AfddNugRY1kpD8xuSvCDJDVV1blU9dow1AQADetKjH553vOSpky5jMHqymQYrOvGxu9+f5P1V9bAkZyS5sqpuSfI/k7y1u78zxhoBgD30pEc/PF8899RJlwFrxorbP6rqoCT/IcnPJflkkj9M8kNJrhxLZf//fU+sqg9X1Rur6sRxvhcAAAxhpT3Z70zy4SQPTvJj3f3s7v7T7v5PSQ5YZrvzq+q2qtq2aHxjVV1fVTdW1aZdvH0nuSvJfkm2r6ReAACYpJVeJ/tN3X3FwoGqemB3/1N3zy2z3QVJXpfkogXbrUvy+iQnZz40b66qy5KsS3LOou1flOTD3X1VVX1fkt9L8pMrrBkAACZipSH7NUmuWDT2scy3iyypu6+uqg2Lho9PcmN335QkVXVJktO6+5wkz1pmd19L8sAV1gsAABOzbMiuqkcmOTTJg6rqiUlq9NJDM986sjsOTXLLguXtSU5YpobnJnlGkgMzf1R8qfXOSnJWkqxfv343SwMAgD23qyPZz8j8yY6HZb5VY4dvJnnFbr5n7WSsl1q5u9+Z5J272ml3n5fkvCSZm5tbcn8AADBuy4bs7r4wyYVV9bzufsdA77k9yeELlg9LcutA+wYAgInbVbvIT3X3W5NsqKpfWfx6d//eTjbblc1Jjq6qI5J8Ocnpmb/RDQAAzIRdXcJv/9HPA5I8ZCePZVXVxZk/QfKYqtpeVWd29z1Jzk7yviTXJbm0u6/dzfoBAGDVqe7Za1+em5vrLVu2TLoMAABmWFVtXepy1rtqF/mj5V7v7l/ck8IAAGAW7erqIlv3ShUAADBDVnJ1EQAA4H7YVbvIH3T3L1XV/8lOrmXd3c8eW2UAADCldtUu8pbRz98ddyEAADArdtUusnX086qq2jfJYzN/RPv67r57L9QHAABTZ1dHspMkVXVqkjcm+Xzmb4t+RFX9fHe/Z5zFAQDANFpRyE7y2iQ/0t03JklVHZXk8iRCNgAALLKrOz7ucNuOgD1yU5LbxlAPAABMvV1dXeS5o6fXVtUVSS7NfE/2v0+yecy1AQDAVNpVu8iPLXj+t0n+zej57UkePpaKAABgyu3q6iI/u7cKAQCAWbHSq4vsl+TMJI9Pst+O8e5+0ZjqAgCAqbXSEx/fkuSRSZ6R5KokhyX55riKAgCAabbSkP2Y7v71JN/q7guTnJrkB8ZXFgAATK+VhuzvjH7eWVXHJXlYkg1jqQgAAKbcSm9Gc15VPTzJrye5LMkBo+cAAMAiKwrZ3f2m0dOrkhw5vnIAAGD6rahdpKoOqqo/rqpPVNXWqvqDqjpo3MUBAMA0WmlP9iWZv43685I8P8kdSf50XEUBAMA0W2lP9iO6+78vWH5NVT1nDPUAAMDUW+mR7A9W1elV9YDR48eTXD7OwgAAYFoteyS7qr6ZpJNUkl9J8tbRSw9IcleSV4+1OgAAmELLhuzufsjeKgQAAGbFSnuyU1XPTvL00eKHuvvd4ykJAACm20ov4Xdukpcm+evR46WjMQAAYJGVHsl+ZpIndPd9SVJVFyb5ZJJN4yoMAACm1UqvLpIkBy54/rCB6wAAgJmx0iPZv5Xkk1X1wcxfaeTpSV4+tqoAAGCK7TJkV9UDktyX5ClJnpz5kP1fu/urY64NAACm0i5DdnffV1Vnd/elSS7bCzUBAMBUW2lP9pVV9Z+r6vCqesSOx1grAwCAKbXSnuwXZf7Oj/9x0fiRw5YDAADTb6Uh+9jMB+ynZT5sfzjJG8dVFAAATLOVhuwLk3wjyR+Nls8Yjf34OIoCAIBpttKQfUx3/+CC5Q9W1afHURAAAEy7lZ74+MmqesqOhao6IclfjqckAACYbis9kn1CkhdW1c2j5fVJrquqzybp7v6XY6kOAACm0EpD9saxVgEAADNkRSG7u7807kIAAGBWrLQnGwAAWCEhGwAABiZkAwDAwIRsAAAY2EqvLsIYHLHp8vSki9iJww7cLx/ZdNJ3l09+7Ydyw+3fmmBFq88hB+ybza88edJlACNvu+bmvOJdn92tbd/xkqfmSY9++MAVjd9zXveRfGr71+/3dg9IctO5pw5f0F6wu5+ZteHAB+2TT736GZMu47scyZ6Q1Rqwk2T7nd/O0879QBIBeym333V3nvyaKyddBpA9C9hJ8rw3fDRbv/S1ASsavz0Jm/clOXLT5cMWtBcI2OzKnf94T57w39436TK+S8iekNUasHf48p3fTpJ8/g4Beym333X3pEsAkrxn21f2eB9/ddPfDVDJ3rPt1m/s0fb3DVTH3rSnn5m14c5/vGfSJXyXkD0hNekCduHQA/dLkhx18P4TrmT1OuSAfSddApDklOMetcf7eMqRBw1Qyd5z3Pc/dI+2n8Z//Pf0M7M2HPig1dMJPY2/ZzPhC+eeumqD9sKe7CtfdmKOPkTQXkxPNqweLzhhfX7r3/3Abm8/jT3Zf3720/KEwx62W9tOa0/2nnxm1obV1pNd3au7caGq/nWSn8z8SZrHdvdTd7XN3Nxcb9myZey1AQCwdlXV1u6e29lrYz2SXVXnV9VtVbVt0fjGqrq+qm6sqk3L7aO7P9zdL07y7iQXjrNeAAAYwrgbVy5I8rokF+0YqKp1SV6f5OQk25NsrqrLkqxLcs6i7V/U3beNnr8gyc+NuV4AANhjYw3Z3X11VW1YNHx8khu7+6YkqapLkpzW3eckedbO9lNV65N8vbudWgwAwKo3iRMfD01yy4Ll7aOx5ZyZ5E+WW6GqzqqqLVW15fbbb9/DEgEAYPdNImTv7KIay5592d2v7u6P7mKd87p7rrvnDjnkkD0qEAAA9sQkQvb2JIcvWD4sya0TqAMAAMZiEiF7c5Kjq+qIqto3yelJLptAHQAAMBbjvoTfxUk+luSYqtpeVWd29z1Jzk7yviTXJbm0u68dZx0AALA3jfvqImcsMX5FkivG+d4AADApbqsOAAADE7IBAGBgQjYAAAxMyAYAgIEJ2QAAMDAhGwAABiZkAwDAwIRsAAAYmJANAAADE7IBAGBgQjYAAAxMyAYAgIEJ2QAAMDAhGwAABiZkAwDAwIRsAAAYmJANAAADE7IBAGBgQjYAAAxMyAYAgIEJ2QAAMDAhGwAABiZkAwDAwIRsAAAYmJANAAADE7IBAGBgQjYAAAxMyAYAgIEJ2QAAMDAhGwAABiZkAwDAwIRsAAAYmJANAAADE7IBAGBgQjYAAAxMyAYAgIEJ2QAAMDAhGwAABiZkAwDAwIRsAAAYmJANAAADE7IBAGBgQjYAAAxMyAYAgIEJ2QAAMDAhGwAABiZkAwDAwIRsAAAY2KoP2VV1bFVdWlVvqKrnT7oeAADYlbGG7Ko6v6puq6pti8Y3VtX1VXVjVW3axW5OSfLH3f2SJC8cW7EAADCQfca8/wuSvC7JRTsGqmpdktcnOTnJ9iSbq+qyJOuSnLNo+xcleUuSV1fVs5McNOZ6AQBgj401ZHf31VW1YdHw8Ulu7O6bkqSqLklyWnefk+RZS+zqF0bh/J1jKxYAAAYy7iPZO3NoklsWLG9PcsJSK49C+iuS7J/kd5ZZ76wkZyXJ+vXrh6gTAAB2yyRCdu1krJdaubu/mFF4Xk53n5fkvCSZm5tbcn8AADBuk7i6yPYkhy9YPizJrROoAwAAxmISIXtzkqOr6oiq2jfJ6Ukum0AdAAAwFuO+hN/FST6W5Jiq2l5VZ3b3PUnOTvK+JNclubS7rx1nHQAAsDeN++oiZywxfkWSK8b53gAAMCmr/o6PAAAwbYRsAAAYmJANAAADE7IBAGBgQjYAAAxMyAYAgIEJ2QAAMDAhGwAABiZkAwDAwIRsAAAYmJANAAADE7IBAGBgQjYAAAxMyAYAgIEJ2QAAMDAhGwAABrbPpAuYFedecV3eePVN93u7ByS56dxThy8IAICJcSR7ALsbsJPkviRHbrp82IIAAJgoIXsA7732q3u0/X0D1QEAwOogZA9g4+MfuUfbmwQAgNki3w1g0zMflxc//cjd2lZPNgDA7KnunnQNg5ubm+stW7ZMugwAAGZYVW3t7rmdveZINgAADEzIBgCAgQnZAAAwMCEbAAAGJmQDAMDAhGwAABiYkA0AAAMTsgEAYGBCNgAADGwm7/hYVbcn+dIE3vrgJHdM4H3Zu8zz2mCe1wbzPPvM8dowqXl+dHcfsrMXZjJkT0pVbVnq1prMDvO8NpjntcE8zz5zvDasxnnWLgIAAAMTsgEAYGBC9rDOm3QB7BXmeW0wz2uDeZ595nhtWHXzrCcbAAAG5kg2AAAMTMgeSFVtrKrrq+rGqto06XpYuao6vKo+WFXXVdW1VfXS0fgjqurKqrph9PPhC7Z5+Wiur6+qZywYf1JVfXb02h9VVU3iM7G0qlpXVZ+sqnePls3zjKmqA6vq7VX1N6Pf6x82z7Olqn559Pf1tqq6uKr2M8fTr6rOr6rbqmrbgrHB5rWqHlhVfzoav6aqNozz8wjZA6iqdUlen+SUJMcmOaOqjp1sVdwP9yR5WXc/LslTkvzCaP42JflAdx+d5AOj5YxeOz3J45NsTPI/Rv8NJMkbkpyV5OjRY+Pe/CCsyEuTXLdg2TzPnj9M8t7ufmySH8z8fJvnGVFVhyb5xSRz3X1cknWZn0NzPP0uyD+fgyHn9cwkX+vuxyT5/SS/PbZPEiF7KMcnubG7b+ruu5NckuS0CdfECnX3V7r7E6Pn38z8P8iHZn4OLxytdmGS54yen5bkku7+p+7+QpIbkxxfVY9K8tDu/ljPn+xw0YJtWAWq6rAkpyZ504Jh8zxDquqhSZ6e5M1J0t13d/edMc+zZp8kD6qqfZI8OMmtMcdTr7uvTvL3i4aHnNeF+3p7kpPG+e2FkD2MQ5PcsmB5+2iMKTP66uiJSa5J8n3d/ZVkPogn+d7RakvN96Gj54vHWT3+IMmvJrlvwZh5ni1HJrk9yZ+M2oLeVFX7xzzPjO7+cpLfTXJzkq8k+Xp3/9+Y41k15Lx+d5vuvifJ15McNK7Chexh7Oz/gly2ZcpU1QFJ3pHkl7r7G8utupOxXmacVaCqnpXktu7eutJNdjJmnle/fZL8UJI3dPcTk3wro6+Xl2Cep8yoJ/e0JEck+f4k+1fVTy23yU7GzPH025153atzLmQPY3uSwxcsH5b5r66YElX1PZkP2P+ru985Gv7b0ddOGf28bTS+1HxvHz1fPM7q8K+SPLuqvpj5lq4fraq3xjzPmu1Jtnf3NaPlt2c+dJvn2fFvk3yhu2/v7u8keWeSp8Ycz6oh5/W724xajR6Wf96eMhghexibkxxdVUdU1b6Zb8S/bMI1sUKjfqw3J7muu39vwUuXJfmZ0fOfSfK/F4yfPjpL+YjMn1Tx8dHXWN+sqqeM9vnCBdswYd398u4+rLs3ZP539C+6+6dinmdKd381yS1Vdcxo6KQkfx3zPEtuTvKUqnrwaG5Oyvy5NOZ4Ng05rwv39fzM/zswvm8vuttjgEeSZyb5XJLPJ/m1Sdfjcb/m7mmZ/7roM0k+NXo8M/N9Wh9IcsPo5yMWbPNro7m+PskpC8bnkmwbvfa6jG745LG6HklOTPLu0XPzPGOPJE9IsmX0O/3nSR5unmfrkeS/Jfmb0fy8JckDzfH0P5JcnPk+++9k/qjzmUPOa5L9kvxZ5k+S/HiSI8f5edzxEQAABqZdBAAABiZkAwDAwIRsAAAYmJANAAADE7IBAGBgQjbAGlRVz6mqYwfa1wVV9fwh9gUwK4RsgLXpOUkGCdkA/HNCNsCUqaoNVfU3VXVhVX2mqt4+uvvdSVX1yar6bFWdX1UPHK1/blX99Wjd362qpyZ5dpLfqapPVdVRVfWYqnp/VX26qj4xGquq+p2q2jba50+M9ldV9brRPi9P8r0LantSVV1VVVur6n07bocMsNa4GQ3AlKmqDUm+kORp3f2XVXV+kpuS/HySk7r7c1V1UZJPJLkoyceSPLa7u6oO7O47q+qCzN/18u2jfV6T5NzufldV7Zf5gzCnJHlxko1JDk6yOckJSX44yUtG49+X+duW/1zmb118VZLTuvv2USh/Rne/aPx/KgCriyPZANPplu7+y9HztyY5KckXuvtzo7ELkzw9yTeSfDvJm6rquUn+YfGOquohSQ7t7nclSXd/u7v/IcnTklzc3fd2999mPkA/ebTfHeO3JvmL0a6OSXJckiur6lNJXpnksIE/N8BU2GfSBQCwW1b0NWR331NVx2c+hJ+e5OwkP7potVpi86XGl3r/SnJtd//wSmoDmGWOZANMp/VVtSPMnpHk/Uk2VNVjRmM/neSqqjogycO6+4okv5TkCaPXv5nkIUnS3d9Isr2qnpMkVfXAqnpwkquT/ERVrauqQzJ/BPvjo/HTR+OPSvIjo31en+SQHXVV1fdU1ePH8ukBVjkhG2A6XZfkZ6rqM0kekeT3k/xskj+rqs8muS/JGzMfpN89Wu+qJL882v6SJP9ldKLkUZkP5b84Wu+jSR6Z5F1JPpPk05lvCfnV7v7qaPyGJJ9N8obRftPddyd5fpLfrqpPJ/lUkqeO8w8BYLVy4iPAlBmd+Pju7j5u0rUAsHOOZAMAwMAcyQYAgIE5kg0AAAMTsgEAYGBCNgAADEzIBgCAgQnZAAAwMCEbAAAG9v8Ay9M1qUiyOwgAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 864x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, axes = plt.subplots(1, figsize=(12, 5))\n",
    "plt.semilogy(samplespace, p4, '.', )\n",
    "axes.set_xlabel('postcode')\n",
    "axes.set_ylabel('probability')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The result is different because the process is different. We are no longer asserting the same constraints as before -- we are only asserting one single constraint. So this will in general have higher entropy (be flatter) than the more constrained model."
   ]
  }
 ],
 "metadata": {
  "@webio": {
   "lastCommId": null,
   "lastKernelId": null
  },
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
